Skip to content

AFNI/NIfTI Server

Sections
Personal tools
You are here: Home » AFNI » Documentation
 

Doxygen Source Code Documentation

Main Page   Alphabetical List   Data Structures   File List   Data Fields   Globals   Search  

poly.c

Go to the documentation of this file.
00001 /*<html><pre>  -<a                             href="qh-poly.htm"
00002   >-------------------------------</a><a name="TOP">-</a>
00003 
00004    poly.c 
00005    implements polygons and simplices
00006 
00007    see qh-poly.htm, poly.h and qhull.h
00008 
00009    infrequent code is in poly2.c 
00010    (all but top 50 and their callers 12/3/95)
00011 
00012    copyright (c) 1993-2001, The Geometry Center
00013 */
00014 
00015 #include "qhull_a.h"
00016 
00017 /*======== functions in alphabetical order ==========*/
00018 
00019 /*-<a                             href="qh-poly.htm#TOC"
00020   >-------------------------------</a><a name="appendfacet">-</a>
00021   
00022   qh_appendfacet( facet )
00023     appends facet to end of qh.facet_list,
00024 
00025   returns:
00026     updates qh.facet_list, facet_tail, newfacet_list, facet_next
00027     increments qh.numfacets
00028   
00029   notes:
00030     assumes qh.facet_list/facet_tail is defined (createsimplex)
00031 
00032 */
00033 void qh_appendfacet(facetT *facet) {
00034   facetT *tail= qh facet_tail;
00035 
00036   if (tail == qh newfacet_list)
00037     qh newfacet_list= facet;
00038   if (tail == qh facet_next)
00039     qh facet_next= facet;
00040   facet->previous= tail->previous;
00041   facet->next= tail;
00042   if (tail->previous)
00043     tail->previous->next= facet;
00044   else
00045     qh facet_list= facet;
00046   tail->previous= facet;
00047   qh num_facets++;
00048   trace4((qh ferr, "qh_appendfacet: append f%d to facet_list\n", facet->id));
00049 } /* appendfacet */
00050 
00051 
00052 /*-<a                             href="qh-poly.htm#TOC"
00053   >-------------------------------</a><a name="appendvertex">-</a>
00054   
00055   qh_appendvertex( vertex )
00056     appends vertex to end of qh.vertex_list,
00057 
00058   returns:
00059     sets vertex->newlist
00060     updates qh.vertex_list, vertex_tail, newvertex_list
00061     increments qh.num_vertices
00062 
00063   notes:
00064     assumes qh.vertex_list/vertex_tail is defined (createsimplex)
00065 
00066 */
00067 void qh_appendvertex (vertexT *vertex) {
00068   vertexT *tail= qh vertex_tail;
00069 
00070   if (tail == qh newvertex_list)
00071     qh newvertex_list= vertex;
00072   vertex->newlist= True;
00073   vertex->previous= tail->previous;
00074   vertex->next= tail;
00075   if (tail->previous)
00076     tail->previous->next= vertex;
00077   else
00078     qh vertex_list= vertex;
00079   tail->previous= vertex;
00080   qh num_vertices++;
00081   trace4((qh ferr, "qh_appendvertex: append v%d to vertex_list\n", vertex->id));
00082 } /* appendvertex */
00083 
00084 
00085 /*-<a                             href="qh-poly.htm#TOC"
00086   >-------------------------------</a><a name="attachnewfacets">-</a>
00087   
00088   qh_attachnewfacets( )
00089     attach horizon facets to new facets in qh.newfacet_list
00090     newfacets have neighbor and ridge links to horizon but not vice versa
00091     only needed for qh.ONLYgood
00092 
00093   returns:
00094     set qh.NEWfacets
00095     horizon facets linked to new facets 
00096       ridges changed from visible facets to new facets
00097       simplicial ridges deleted
00098     qh.visible_list, no ridges valid
00099     facet->f.replace is a newfacet (if any)
00100 
00101   design:
00102     delete interior ridges and neighbor sets by
00103       for each visible, non-simplicial facet
00104         for each ridge
00105           if last visit or if neighbor is simplicial
00106             if horizon neighbor
00107               delete ridge for horizon's ridge set
00108             delete ridge
00109         erase neighbor set
00110     attach horizon facets and new facets by
00111       for all new facets
00112         if corresponding horizon facet is simplicial
00113           locate corresponding visible facet {may be more than one}
00114           link visible facet to new facet
00115           replace visible facet with new facet in horizon
00116         else it's non-simplicial
00117           for all visible neighbors of the horizon facet
00118             link visible neighbor to new facet
00119             delete visible neighbor from horizon facet
00120           append new facet to horizon's neighbors
00121           the first ridge of the new facet is the horizon ridge
00122           link the new facet into the horizon ridge
00123 */
00124 void qh_attachnewfacets (void ) {
00125   facetT *newfacet= NULL, *neighbor, **neighborp, *horizon, *visible;
00126   ridgeT *ridge, **ridgep;
00127 
00128   qh NEWfacets= True;
00129   trace3((qh ferr, "qh_attachnewfacets: delete interior ridges\n"));
00130   qh visit_id++;
00131   FORALLvisible_facets {
00132     visible->visitid= qh visit_id;
00133     if (visible->ridges) {
00134       FOREACHridge_(visible->ridges) {
00135         neighbor= otherfacet_(ridge, visible);
00136         if (neighbor->visitid == qh visit_id
00137             || (!neighbor->visible && neighbor->simplicial)) {
00138           if (!neighbor->visible)  /* delete ridge for simplicial horizon */
00139             qh_setdel (neighbor->ridges, ridge);
00140           qh_setfree (&(ridge->vertices)); /* delete on 2nd visit */
00141           qh_memfree (ridge, sizeof(ridgeT));
00142         }
00143       }
00144       SETfirst_(visible->ridges)= NULL;
00145     }
00146     SETfirst_(visible->neighbors)= NULL;
00147   }
00148   trace1((qh ferr, "qh_attachnewfacets: attach horizon facets to new facets\n"));
00149   FORALLnew_facets {
00150     horizon= SETfirstt_(newfacet->neighbors, facetT);
00151     if (horizon->simplicial) {
00152       visible= NULL;
00153       FOREACHneighbor_(horizon) {   /* may have more than one horizon ridge */
00154         if (neighbor->visible) {
00155           if (visible) {
00156             if (qh_setequal_skip (newfacet->vertices, 0, horizon->vertices,
00157                                   SETindex_(horizon->neighbors, neighbor))) {
00158               visible= neighbor;
00159               break;
00160             }
00161           }else
00162             visible= neighbor;
00163         }
00164       }
00165       if (visible) {
00166         visible->f.replace= newfacet;
00167         qh_setreplace (horizon->neighbors, visible, newfacet);
00168       }else {
00169         fprintf (qh ferr, "qhull internal error (qh_attachnewfacets): couldn't find visible facet for horizon f%d of newfacet f%d\n",
00170                  horizon->id, newfacet->id);
00171         qh_errexit2 (qh_ERRqhull, horizon, newfacet);
00172       }
00173     }else { /* non-simplicial, with a ridge for newfacet */
00174       FOREACHneighbor_(horizon) {    /* may hold for many new facets */
00175         if (neighbor->visible) {
00176           neighbor->f.replace= newfacet;
00177           qh_setdelnth (horizon->neighbors,
00178                         SETindex_(horizon->neighbors, neighbor));
00179           neighborp--; /* repeat */
00180         }
00181       }
00182       qh_setappend (&horizon->neighbors, newfacet);
00183       ridge= SETfirstt_(newfacet->ridges, ridgeT);
00184       if (ridge->top == horizon)
00185         ridge->bottom= newfacet;
00186       else
00187         ridge->top= newfacet;
00188       }
00189   } /* newfacets */
00190   if (qh PRINTstatistics) {
00191     FORALLvisible_facets {
00192       if (!visible->f.replace) 
00193         zinc_(Zinsidevisible);
00194     }
00195   }
00196 } /* attachnewfacets */
00197 
00198 /*-<a                             href="qh-poly.htm#TOC"
00199   >-------------------------------</a><a name="checkflipped">-</a>
00200   
00201   qh_checkflipped( facet, dist, allerror )
00202     checks facet orientation to interior point
00203 
00204     if allerror set,
00205       tests against qh.DISTround
00206     else
00207       tests against 0 since tested against DISTround before
00208 
00209   returns:
00210     False if it flipped orientation (sets facet->flipped)
00211     distance if non-NULL
00212 */
00213 boolT qh_checkflipped (facetT *facet, realT *distp, boolT allerror) {
00214   realT dist;
00215 
00216   if (facet->flipped && !distp)
00217     return False;
00218   zzinc_(Zdistcheck);
00219   qh_distplane(qh interior_point, facet, &dist);
00220   if (distp)
00221     *distp= dist;
00222   if ((allerror && dist > -qh DISTround)|| (!allerror && dist >= 0.0)) {
00223     facet->flipped= True;
00224     zzinc_(Zflippedfacets);
00225     trace0((qh ferr, "qh_checkflipped: facet f%d is flipped, distance= %6.12g during p%d\n",
00226               facet->id, dist, qh furthest_id));
00227     qh_precision ("flipped facet");
00228     return False;
00229   }
00230   return True;
00231 } /* checkflipped */
00232 
00233 /*-<a                             href="qh-poly.htm#TOC"
00234   >-------------------------------</a><a name="delfacet">-</a>
00235   
00236   qh_delfacet( facet )
00237     removes facet from facet_list and frees up its memory
00238 
00239   notes:
00240     assumes vertices and ridges already freed
00241 */
00242 void qh_delfacet(facetT *facet) {
00243   void **freelistp; /* used !qh_NOmem */
00244 
00245   trace5((qh ferr, "qh_delfacet: delete f%d\n", facet->id));
00246   if (facet == qh tracefacet)
00247     qh tracefacet= NULL;
00248   if (facet == qh GOODclosest)
00249     qh GOODclosest= NULL;
00250   qh_removefacet(facet);
00251   qh_memfree_(facet->normal, qh normal_size, freelistp);
00252   if (qh CENTERtype == qh_ASvoronoi) {   /* uses macro calls */
00253     qh_memfree_(facet->center, qh center_size, freelistp);
00254   }else /* AScentrum */ {
00255     qh_memfree_(facet->center, qh normal_size, freelistp);
00256   }
00257   qh_setfree(&(facet->neighbors));
00258   if (facet->ridges)
00259     qh_setfree(&(facet->ridges));
00260   qh_setfree(&(facet->vertices));
00261   if (facet->outsideset)
00262     qh_setfree(&(facet->outsideset));
00263   if (facet->coplanarset)
00264     qh_setfree(&(facet->coplanarset));
00265   qh_memfree_(facet, sizeof(facetT), freelistp);
00266 } /* delfacet */
00267 
00268 
00269 /*-<a                             href="qh-poly.htm#TOC"
00270   >-------------------------------</a><a name="deletevisible">-</a>
00271   
00272   qh_deletevisible()
00273     delete visible facets and vertices
00274 
00275   returns:
00276     deletes each facet and removes from facetlist
00277     at exit, qh.visible_list empty (== qh.newfacet_list)
00278 
00279   notes:
00280     ridges already deleted
00281     horizon facets do not reference facets on qh.visible_list
00282     new facets in qh.newfacet_list
00283     uses   qh.visit_id;
00284 */
00285 void qh_deletevisible (void /*qh visible_list*/) {
00286   facetT *visible, *nextfacet;
00287   vertexT *vertex, **vertexp;
00288   int numvisible= 0, numdel= qh_setsize(qh del_vertices);
00289 
00290   trace1((qh ferr, "qh_deletevisible: delete %d visible facets and %d vertices\n",
00291          qh num_visible, numdel));
00292   for (visible= qh visible_list; visible && visible->visible; 
00293                 visible= nextfacet) { /* deleting current */
00294     nextfacet= visible->next;        
00295     numvisible++;
00296     qh_delfacet(visible);
00297   }
00298   if (numvisible != qh num_visible) {
00299     fprintf (qh ferr, "qhull internal error (qh_deletevisible): qh num_visible %d is not number of visible facets %d\n",
00300              qh num_visible, numvisible);
00301     qh_errexit (qh_ERRqhull, NULL, NULL);
00302   }
00303   qh num_visible= 0;
00304   zadd_(Zvisfacettot, numvisible);
00305   zmax_(Zvisfacetmax, numvisible);
00306   zadd_(Zdelvertextot, numdel);
00307   zmax_(Zdelvertexmax, numdel);
00308   FOREACHvertex_(qh del_vertices) 
00309     qh_delvertex (vertex);
00310   qh_settruncate (qh del_vertices, 0);
00311 } /* deletevisible */
00312 
00313 /*-<a                             href="qh-poly.htm#TOC"
00314   >-------------------------------</a><a name="facetintersect">-</a>
00315   
00316   qh_facetintersect( facetA, facetB, skipa, skipB, prepend )
00317     return vertices for intersection of two simplicial facets
00318     may include 1 prepended entry (if more, need to settemppush)
00319     
00320   returns:
00321     returns set of qh.hull_dim-1 + prepend vertices
00322     returns skipped index for each test and checks for exactly one
00323 
00324   notes:
00325     does not need settemp since set in quick memory
00326   
00327   see also:
00328     qh_vertexintersect and qh_vertexintersect_new
00329     use qh_setnew_delnthsorted to get nth ridge (no skip information)
00330 
00331   design:
00332     locate skipped vertex by scanning facet A's neighbors
00333     locate skipped vertex by scanning facet B's neighbors
00334     intersect the vertex sets
00335 */
00336 setT *qh_facetintersect (facetT *facetA, facetT *facetB,
00337                          int *skipA,int *skipB, int prepend) {
00338   setT *intersect;
00339   int dim= qh hull_dim, i, j;
00340   facetT **neighborsA, **neighborsB;
00341 
00342   neighborsA= SETaddr_(facetA->neighbors, facetT);
00343   neighborsB= SETaddr_(facetB->neighbors, facetT);
00344   i= j= 0;
00345   if (facetB == *neighborsA++)
00346     *skipA= 0;
00347   else if (facetB == *neighborsA++)
00348     *skipA= 1;
00349   else if (facetB == *neighborsA++)
00350     *skipA= 2;
00351   else {
00352     for (i= 3; i < dim; i++) {
00353       if (facetB == *neighborsA++) {
00354         *skipA= i;
00355         break;
00356       }
00357     }
00358   }
00359   if (facetA == *neighborsB++)
00360     *skipB= 0;
00361   else if (facetA == *neighborsB++)
00362     *skipB= 1;
00363   else if (facetA == *neighborsB++)
00364     *skipB= 2;
00365   else {
00366     for (j= 3; j < dim; j++) {
00367       if (facetA == *neighborsB++) {
00368         *skipB= j;
00369         break;
00370       }
00371     }
00372   }
00373   if (i >= dim || j >= dim) {
00374     fprintf (qh ferr, "qhull internal error (qh_facetintersect): f%d or f%d not in others neighbors\n",
00375             facetA->id, facetB->id);
00376     qh_errexit2 (qh_ERRqhull, facetA, facetB);
00377   }
00378   intersect= qh_setnew_delnthsorted (facetA->vertices, qh hull_dim, *skipA, prepend);
00379   trace4((qh ferr, "qh_facetintersect: f%d skip %d matches f%d skip %d\n",
00380           facetA->id, *skipA, facetB->id, *skipB));
00381   return(intersect);
00382 } /* facetintersect */
00383 
00384 /*-<a                             href="qh-poly.htm#TOC"
00385   >-------------------------------</a><a name="gethash">-</a>
00386   
00387   qh_gethash( hashsize, set, size, firstindex, skipelem )
00388     return hashvalue for a set with firstindex and skipelem
00389 
00390   notes:
00391     assumes at least firstindex+1 elements
00392     assumes skipelem is NULL, in set, or part of hash
00393     
00394     hashes memory addresses which may change over different runs of the same data
00395     using sum for hash does badly in high d
00396 */
00397 unsigned qh_gethash (int hashsize, setT *set, int size, int firstindex, void *skipelem) {
00398   void **elemp= SETelemaddr_(set, firstindex, void);
00399   ptr_intT hash = 0, elem;
00400   int i;
00401 
00402   switch (size-firstindex) {
00403   case 1:
00404     hash= (ptr_intT)(*elemp) - (ptr_intT) skipelem;
00405     break;
00406   case 2:
00407     hash= (ptr_intT)(*elemp) + (ptr_intT)elemp[1] - (ptr_intT) skipelem;
00408     break;
00409   case 3:
00410     hash= (ptr_intT)(*elemp) + (ptr_intT)elemp[1] + (ptr_intT)elemp[2]
00411       - (ptr_intT) skipelem;
00412     break;
00413   case 4:
00414     hash= (ptr_intT)(*elemp) + (ptr_intT)elemp[1] + (ptr_intT)elemp[2]
00415       + (ptr_intT)elemp[3] - (ptr_intT) skipelem;
00416     break;
00417   case 5:
00418     hash= (ptr_intT)(*elemp) + (ptr_intT)elemp[1] + (ptr_intT)elemp[2]
00419       + (ptr_intT)elemp[3] + (ptr_intT)elemp[4] - (ptr_intT) skipelem;
00420     break;
00421   case 6:
00422     hash= (ptr_intT)(*elemp) + (ptr_intT)elemp[1] + (ptr_intT)elemp[2]
00423       + (ptr_intT)elemp[3] + (ptr_intT)elemp[4]+ (ptr_intT)elemp[5]
00424       - (ptr_intT) skipelem;
00425     break;
00426   default:
00427     hash= 0;
00428     i= 3;
00429     do {     /* this is about 10% in 10-d */
00430       if ((elem= (ptr_intT)*elemp++) != (ptr_intT)skipelem) {
00431         hash ^= (elem << i) + (elem >> (32-i));
00432         i += 3;
00433         if (i >= 32)
00434           i -= 32;
00435       }
00436     }while(*elemp);
00437     break;
00438   }
00439   hash %= (ptr_intT) hashsize;
00440   /* hash= 0; for debugging purposes */
00441   return hash;
00442 } /* gethash */
00443 
00444 /*-<a                             href="qh-poly.htm#TOC"
00445   >-------------------------------</a><a name="makenewfacet">-</a>
00446   
00447   qh_makenewfacet( vertices, toporient, horizon )
00448     creates a toporient? facet from vertices
00449 
00450   returns:
00451     returns newfacet
00452       adds newfacet to qh.facet_list 
00453       newfacet->neighbor= horizon, but not vice versa
00454       newfacet->vertices= vertices
00455     newvertex_list updated with vertices
00456 */
00457 facetT *qh_makenewfacet(setT *vertices, boolT toporient,facetT *horizon) {
00458   facetT *newfacet;
00459   vertexT *vertex, **vertexp;
00460 
00461   FOREACHvertex_(vertices) {
00462     if (!vertex->newlist) {
00463       qh_removevertex (vertex);
00464       qh_appendvertex (vertex);
00465     }
00466   }
00467   newfacet= qh_newfacet();
00468   newfacet->vertices= vertices;
00469   newfacet->toporient= toporient;
00470   qh_setappend(&(newfacet->neighbors), horizon);
00471   qh_appendfacet(newfacet);
00472   return(newfacet);
00473 } /* makenewfacet */
00474 
00475 
00476 /*-<a                             href="qh-poly.htm#TOC"
00477   >-------------------------------</a><a name="makenewplanes">-</a>
00478   
00479   qh_makenewplanes()
00480     make new hyperplanes for facets on qh.newfacet_list
00481 
00482   returns:
00483     all facets have hyperplanes or are marked for   merging
00484     doesn't create hyperplane if horizon is coplanar (will merge)
00485     updates qh.min_vertex if qh.JOGGLEmax
00486 
00487   notes:
00488     facet->f.samecycle is defined for facet->mergehorizon facets
00489 */
00490 void qh_makenewplanes (void /* newfacet_list */) {
00491   facetT *newfacet;
00492 
00493   FORALLnew_facets {
00494     if (!newfacet->mergehorizon)
00495       qh_setfacetplane (newfacet);  
00496   }
00497   if (qh JOGGLEmax < REALmax/2)  
00498     minimize_(qh min_vertex, -wwval_(Wnewvertexmax));
00499 } /* makenewplanes */
00500 
00501 /*-<a                             href="qh-poly.htm#TOC"
00502   >-------------------------------</a><a name="makenew_nonsimplicial">-</a>
00503   
00504   qh_makenew_nonsimplicial( visible, apex, numnew )
00505     make new facets for ridges of a visible facet
00506     
00507   returns:
00508     first newfacet, bumps numnew as needed
00509     attaches new facets if !qh.ONLYgood
00510     marks ridge neighbors for simplicial visible
00511     if (qh.ONLYgood)
00512       ridges on newfacet, horizon, and visible
00513     else
00514       ridge and neighbors between newfacet and   horizon
00515       visible facet's ridges are deleted    
00516 
00517   notes:
00518     qh.visit_id if visible has already been processed
00519     sets neighbor->seen for building f.samecycle
00520       assumes all 'seen' flags initially false
00521     
00522   design:
00523     for each ridge of visible facet
00524       get neighbor of visible facet
00525       if neighbor was already processed
00526         delete the ridge (will delete all visible facets later)
00527       if neighbor is a horizon facet
00528         create a new facet
00529         if neighbor coplanar
00530           adds newfacet to f.samecycle for later merging
00531         else 
00532           updates neighbor's neighbor set
00533           (checks for non-simplicial facet with multiple ridges to visible facet)
00534         updates neighbor's ridge set
00535         (checks for simplicial neighbor to non-simplicial visible facet)
00536           
00537 */
00538 #ifndef qh_NOmerge
00539 facetT *qh_makenew_nonsimplicial (facetT *visible, vertexT *apex, int *numnew) {
00540   void **freelistp; /* used !qh_NOmem */
00541   ridgeT *ridge, **ridgep;
00542   facetT *neighbor, *newfacet= NULL, *samecycle;
00543   setT *vertices;
00544   boolT toporient;
00545   int ridgeid;
00546 
00547   FOREACHridge_(visible->ridges) {
00548     ridgeid= ridge->id;
00549     neighbor= otherfacet_(ridge, visible);
00550     if (neighbor->visible) {
00551       if (!qh ONLYgood) {
00552         if (neighbor->visitid == qh visit_id) {
00553           qh_setfree (&(ridge->vertices));  /* delete on 2nd visit */
00554           qh_memfree_(ridge, sizeof(ridgeT), freelistp);
00555         }
00556       }
00557     }else {  /* neighbor is an horizon facet */
00558       toporient= (ridge->top == visible);
00559       vertices= qh_setnew (qh hull_dim); /* makes sure this is quick */
00560       qh_setappend (&vertices, apex);
00561       qh_setappend_set (&vertices, ridge->vertices);
00562       newfacet= qh_makenewfacet(vertices, toporient, neighbor);
00563       (*numnew)++;
00564       if (neighbor->coplanar) {
00565         newfacet->mergehorizon= True;
00566         if (!neighbor->seen) {
00567           newfacet->f.samecycle= newfacet;
00568           neighbor->f.newcycle= newfacet;
00569         }else {
00570           samecycle= neighbor->f.newcycle;
00571           newfacet->f.samecycle= samecycle->f.samecycle;
00572           samecycle->f.samecycle= newfacet;
00573         }
00574       }
00575       if (qh ONLYgood) {
00576         if (!neighbor->simplicial)
00577           qh_setappend(&(newfacet->ridges), ridge);
00578       }else {  /* qh_attachnewfacets */
00579         if (neighbor->seen) {
00580           if (neighbor->simplicial) {
00581             fprintf (qh ferr, "qhull internal error (qh_makenew_nonsimplicial): simplicial f%d sharing two ridges with f%d\n", 
00582                    neighbor->id, visible->id);
00583             qh_errexit2 (qh_ERRqhull, neighbor, visible);
00584           }
00585           qh_setappend (&(neighbor->neighbors), newfacet);
00586         }else
00587           qh_setreplace (neighbor->neighbors, visible, newfacet);
00588         if (neighbor->simplicial) {
00589           qh_setdel (neighbor->ridges, ridge);
00590           qh_setfree (&(ridge->vertices)); 
00591           qh_memfree (ridge, sizeof(ridgeT));
00592         }else {
00593           qh_setappend(&(newfacet->ridges), ridge);
00594           if (toporient)
00595             ridge->top= newfacet;
00596           else
00597             ridge->bottom= newfacet;
00598         }
00599       trace4((qh ferr, "qh_makenew_nonsimplicial: created facet f%d from v%d and r%d of horizon f%d\n",
00600             newfacet->id, apex->id, ridgeid, neighbor->id));
00601       }
00602     }
00603     neighbor->seen= True;        
00604   } /* for each ridge */
00605   if (!qh ONLYgood)
00606     SETfirst_(visible->ridges)= NULL;
00607   return newfacet;
00608 } /* makenew_nonsimplicial */
00609 #else /* qh_NOmerge */
00610 facetT *qh_makenew_nonsimplicial (facetT *visible, vertexT *apex, int *numnew) {
00611   return NULL;
00612 }
00613 #endif /* qh_NOmerge */
00614 
00615 /*-<a                             href="qh-poly.htm#TOC"
00616   >-------------------------------</a><a name="makenew_simplicial">-</a>
00617   
00618   qh_makenew_simplicial( visible, apex, numnew )
00619     make new facets for simplicial visible facet and apex
00620 
00621   returns:
00622     attaches new facets if (!qh.ONLYgood)
00623       neighbors between newfacet and horizon
00624 
00625   notes:
00626     nop if neighbor->seen or neighbor->visible (see qh_makenew_nonsimplicial)
00627 
00628   design:
00629     locate neighboring horizon facet for visible facet
00630     determine vertices and orientation
00631     create new facet
00632     if coplanar,
00633       add new facet to f.samecycle
00634     update horizon facet's neighbor list        
00635 */
00636 facetT *qh_makenew_simplicial (facetT *visible, vertexT *apex, int *numnew) {
00637   facetT *neighbor, **neighborp, *newfacet= NULL;
00638   setT *vertices;
00639   boolT flip, toporient;
00640   int horizonskip, visibleskip;
00641 
00642   FOREACHneighbor_(visible) {
00643     if (!neighbor->seen && !neighbor->visible) {
00644       vertices= qh_facetintersect(neighbor,visible, &horizonskip, &visibleskip, 1);
00645       SETfirst_(vertices)= apex;
00646       flip= ((horizonskip & 0x1) ^ (visibleskip & 0x1));
00647       if (neighbor->toporient)         
00648         toporient= horizonskip & 0x1;
00649       else
00650         toporient= (horizonskip & 0x1) ^ 0x1;
00651       newfacet= qh_makenewfacet(vertices, toporient, neighbor);
00652       (*numnew)++;
00653       if (neighbor->coplanar && (qh PREmerge || qh MERGEexact)) {
00654 #ifndef qh_NOmerge
00655         newfacet->f.samecycle= newfacet;
00656         newfacet->mergehorizon= True;
00657 #endif
00658       }
00659       if (!qh ONLYgood)
00660         SETelem_(neighbor->neighbors, horizonskip)= newfacet;
00661       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",
00662             newfacet->id, toporient, apex->id, neighbor->id, horizonskip,
00663               neighbor->toporient, visible->id, visibleskip, flip));
00664     }
00665   }
00666   return newfacet;
00667 } /* makenew_simplicial */
00668 
00669 /*-<a                             href="qh-poly.htm#TOC"
00670   >-------------------------------</a><a name="matchneighbor">-</a>
00671   
00672   qh_matchneighbor( newfacet, newskip, hashsize, hashcount )
00673     either match subridge of newfacet with neighbor or add to hash_table
00674 
00675   returns:
00676     duplicate ridges are unmatched and marked by qh_DUPLICATEridge
00677 
00678   notes:
00679     ridge is newfacet->vertices w/o newskip vertex
00680     do not allocate memory (need to free hash_table cleanly)
00681     uses linear hash chains
00682   
00683   see also:
00684     qh_matchduplicates
00685 
00686   design:
00687     for each possible matching facet in qh.hash_table
00688       if vertices match
00689         set ismatch, if facets have opposite orientation
00690         if ismatch and matching facet doesn't have a match
00691           match the facets by updating their neighbor sets
00692         else
00693           indicate a duplicate ridge
00694           set facet hyperplane for later testing
00695           add facet to hashtable
00696           unless the other facet was already a duplicate ridge
00697             mark both facets with a duplicate ridge
00698             add other facet (if defined) to hash table
00699 */
00700 void qh_matchneighbor (facetT *newfacet, int newskip, int hashsize, int *hashcount) {
00701   boolT newfound= False;   /* True, if new facet is already in hash chain */
00702   boolT same, ismatch;
00703   int hash, scan;
00704   facetT *facet, *matchfacet;
00705   int skip, matchskip;
00706 
00707   hash= (int)qh_gethash (hashsize, newfacet->vertices, qh hull_dim, 1, 
00708                      SETelem_(newfacet->vertices, newskip));
00709   trace4((qh ferr, "qh_matchneighbor: newfacet f%d skip %d hash %d hashcount %d\n",
00710           newfacet->id, newskip, hash, *hashcount));
00711   zinc_(Zhashlookup);
00712   for (scan= hash; (facet= SETelemt_(qh hash_table, scan, facetT)); 
00713        scan= (++scan >= hashsize ? 0 : scan)) {
00714     if (facet == newfacet) {
00715       newfound= True;
00716       continue;
00717     }
00718     zinc_(Zhashtests);
00719     if (qh_matchvertices (1, newfacet->vertices, newskip, facet->vertices, &skip, &same)) {
00720       if (SETelem_(newfacet->vertices, newskip) == 
00721           SETelem_(facet->vertices, skip)) {
00722         qh_precision ("two facets with the same vertices");
00723         fprintf (qh ferr, "qhull precision error: Vertex sets are the same for f%d and f%d.  Can not force output.\n",
00724           facet->id, newfacet->id);
00725         qh_errexit2 (qh_ERRprec, facet, newfacet);
00726       }
00727       ismatch= (same == (newfacet->toporient ^ facet->toporient));
00728       matchfacet= SETelemt_(facet->neighbors, skip, facetT);
00729       if (ismatch && !matchfacet) {
00730         SETelem_(facet->neighbors, skip)= newfacet;
00731         SETelem_(newfacet->neighbors, newskip)= facet;
00732         (*hashcount)--;
00733         trace4((qh ferr, "qh_matchneighbor: f%d skip %d matched with new f%d skip %d\n",
00734            facet->id, skip, newfacet->id, newskip));
00735         return;
00736       }
00737       if (!qh PREmerge && !qh MERGEexact) {
00738         qh_precision ("a ridge with more than two neighbors");
00739         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",
00740                  facet->id, newfacet->id, getid_(matchfacet));
00741         qh_errexit2 (qh_ERRprec, facet, newfacet);
00742       }
00743       SETelem_(newfacet->neighbors, newskip)= qh_DUPLICATEridge;
00744       newfacet->dupridge= True;
00745       if (!newfacet->normal)
00746         qh_setfacetplane (newfacet);
00747       qh_addhash (newfacet, qh hash_table, hashsize, hash);
00748       (*hashcount)++;
00749       if (!facet->normal)
00750         qh_setfacetplane (facet);
00751       if (matchfacet != qh_DUPLICATEridge) {
00752         SETelem_(facet->neighbors, skip)= qh_DUPLICATEridge;
00753         facet->dupridge= True;
00754         if (!facet->normal)
00755           qh_setfacetplane (facet);
00756         if (matchfacet) {
00757           matchskip= qh_setindex (matchfacet->neighbors, facet);
00758           SETelem_(matchfacet->neighbors, matchskip)= qh_DUPLICATEridge;
00759           matchfacet->dupridge= True;
00760           if (!matchfacet->normal)
00761             qh_setfacetplane (matchfacet);
00762           qh_addhash (matchfacet, qh hash_table, hashsize, hash);
00763           *hashcount += 2;
00764         }
00765       }
00766       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",
00767            newfacet->id, newskip, facet->id, skip, 
00768            (matchfacet == qh_DUPLICATEridge ? -2 : getid_(matchfacet)), 
00769            ismatch, hash));
00770       return; /* end of duplicate ridge */
00771     }
00772   }
00773   if (!newfound) 
00774     SETelem_(qh hash_table, scan)= newfacet;  /* same as qh_addhash */
00775   (*hashcount)++;
00776   trace4((qh ferr, "qh_matchneighbor: no match for f%d skip %d at hash %d\n",
00777            newfacet->id, newskip, hash));
00778 } /* matchneighbor */
00779 
00780 
00781 /*-<a                             href="qh-poly.htm#TOC"
00782   >-------------------------------</a><a name="matchnewfacets">-</a>
00783   
00784   qh_matchnewfacets()
00785     match newfacets in qh.newfacet_list to their newfacet neighbors
00786 
00787   returns:
00788     qh.newfacet_list with full neighbor sets
00789       get vertices with nth neighbor by deleting nth vertex
00790     if qh.PREmerge/MERGEexact or qh.FORCEoutput 
00791       all facets check for flipped (also prevents point partitioning)
00792     if duplicate ridges and qh.PREmerge/MERGEexact
00793       facet->dupridge set
00794       missing neighbor links identifies extra ridges to be merging
00795 
00796   notes:
00797     newfacets already have neighbor[0] (horizon facet)
00798     assumes qh.hash_table is NULL
00799     vertex->neighbors has not been updated yet
00800     do not allocate memory after qh.hash_table (need to free it cleanly)
00801 
00802   design:
00803     delete neighbor sets for all new facets
00804     initialize a hash table
00805     for all new facets
00806       match facet with neighbors
00807     if unmatched facets (due to duplicate ridges)
00808       for each new facet with a duplicate ridge
00809         match it with a facet
00810     check for flipped facets
00811 */
00812 void qh_matchnewfacets (void) {
00813   int numnew=0, hashcount=0, newskip;
00814   facetT *newfacet, *neighbor;
00815   int dim= qh hull_dim, hashsize, neighbor_i, neighbor_n;
00816   setT *neighbors;
00817 #ifndef qh_NOtrace
00818   int facet_i, facet_n, numfree= 0;
00819   facetT *facet;
00820 #endif
00821   
00822   trace1((qh ferr, "qh_matchnewfacets: match neighbors for new facets.\n"));
00823   FORALLnew_facets {
00824     numnew++;
00825     {  /* inline qh_setzero (newfacet->neighbors, 1, qh hull_dim); */
00826       neighbors= newfacet->neighbors;
00827       neighbors->e[neighbors->maxsize].i= dim+1; /*may be overwritten*/
00828       memset ((char *)SETelemaddr_(neighbors, 1, void), 0, dim * SETelemsize);
00829     }    
00830   }
00831   qh_newhashtable (numnew*(qh hull_dim-1)); /* twice what is normally needed,
00832                                      but every ridge could be DUPLICATEridge */
00833   hashsize= qh_setsize (qh hash_table);
00834   FORALLnew_facets {
00835     for (newskip=1; newskip<qh hull_dim; newskip++) /* furthest/horizon already matched */
00836       qh_matchneighbor (newfacet, newskip, hashsize, &hashcount);
00837 #if 0   /* use the following to trap hashcount errors */
00838     {
00839       int count= 0, k;
00840       facetT *facet, *neighbor;
00841 
00842       count= 0;
00843       FORALLfacet_(qh newfacet_list) {  /* newfacet already in use */
00844         for (k=1; k < qh hull_dim; k++) {
00845           neighbor= SETelemt_(facet->neighbors, k, facetT);
00846           if (!neighbor || neighbor == qh_DUPLICATEridge)
00847             count++;
00848         }
00849         if (facet == newfacet)
00850           break;
00851       }
00852       if (count != hashcount) {
00853         fprintf (qh ferr, "qh_matchnewfacets: after adding facet %d, hashcount %d != count %d\n",
00854                  newfacet->id, hashcount, count);
00855         qh_errexit (qh_ERRqhull, newfacet, NULL);
00856       }
00857     }
00858 #endif  /* end of trap code */
00859   }
00860   if (hashcount) {
00861     FORALLnew_facets {
00862       if (newfacet->dupridge) {
00863         FOREACHneighbor_i_(newfacet) {
00864           if (neighbor == qh_DUPLICATEridge) {
00865             qh_matchduplicates (newfacet, neighbor_i, hashsize, &hashcount);
00866                     /* this may report MERGEfacet */
00867           }
00868         }
00869       }
00870     }
00871   }
00872   if (hashcount) {
00873     fprintf (qh ferr, "qhull internal error (qh_matchnewfacets): %d neighbors did not match up\n",
00874         hashcount);
00875     qh_printhashtable (qh ferr);
00876     qh_errexit (qh_ERRqhull, NULL, NULL);
00877   }
00878 #ifndef qh_NOtrace
00879   if (qh IStracing >= 2) {
00880     FOREACHfacet_i_(qh hash_table) {
00881       if (!facet)
00882         numfree++;
00883     }
00884     fprintf (qh ferr, "qh_matchnewfacets: %d new facets, %d unused hash entries .  hashsize %d\n",
00885              numnew, numfree, qh_setsize (qh hash_table));
00886   }
00887 #endif /* !qh_NOtrace */
00888   qh_setfree (&qh hash_table);
00889   if (qh PREmerge || qh MERGEexact) {
00890     if (qh IStracing >= 4)
00891       qh_printfacetlist (qh newfacet_list, NULL, qh_ALL);
00892     FORALLnew_facets {
00893       if (newfacet->normal)
00894         qh_checkflipped (newfacet, NULL, qh_ALL);
00895     }
00896   }else if (qh FORCEoutput)
00897     qh_checkflipped_all (qh newfacet_list);  /* prints warnings for flipped */
00898 } /* matchnewfacets */
00899 
00900     
00901 /*-<a                             href="qh-poly.htm#TOC"
00902   >-------------------------------</a><a name="matchvertices">-</a>
00903   
00904   qh_matchvertices( firstindex, verticesA, skipA, verticesB, skipB, same )
00905     tests whether vertices match with a single skip
00906     starts match at firstindex since all new facets have a common vertex
00907 
00908   returns:
00909     true if matched vertices
00910     skip index for each set
00911     sets same iff vertices have the same orientation
00912 
00913   notes:
00914     assumes skipA is in A and both sets are the same size
00915 
00916   design:
00917     set up pointers
00918     scan both sets checking for a match
00919     test orientation
00920 */
00921 boolT qh_matchvertices (int firstindex, setT *verticesA, int skipA, 
00922        setT *verticesB, int *skipB, boolT *same) {
00923   vertexT **elemAp, **elemBp, **skipBp=NULL, **skipAp;
00924 
00925   elemAp= SETelemaddr_(verticesA, firstindex, vertexT);
00926   elemBp= SETelemaddr_(verticesB, firstindex, vertexT);
00927   skipAp= SETelemaddr_(verticesA, skipA, vertexT);
00928   do if (elemAp != skipAp) {
00929     while (*elemAp != *elemBp++) {
00930       if (skipBp)
00931         return False;
00932       skipBp= elemBp;  /* one extra like FOREACH */
00933     }
00934   }while(*(++elemAp));
00935   if (!skipBp)
00936     skipBp= ++elemBp;
00937   *skipB= SETindex_(verticesB, skipB);
00938   *same= !(((ptr_intT)skipA & 0x1) ^ ((ptr_intT)*skipB & 0x1));
00939   trace4((qh ferr, "qh_matchvertices: matched by skip %d (v%d) and skip %d (v%d) same? %d\n",
00940           skipA, (*skipAp)->id, *skipB, (*(skipBp-1))->id, *same));
00941   return (True);
00942 } /* matchvertices */
00943 
00944 /*-<a                             href="qh-poly.htm#TOC"
00945   >-------------------------------</a><a name="newfacet">-</a>
00946   
00947   qh_newfacet()
00948     return a new facet 
00949 
00950   returns:
00951     all fields initialized or cleared   (NULL)
00952     preallocates neighbors set
00953 */
00954 facetT *qh_newfacet(void) {
00955   facetT *facet;
00956   void **freelistp; /* used !qh_NOmem */
00957   
00958   qh_memalloc_(sizeof(facetT), freelistp, facet, facetT);
00959   memset ((char *)facet, 0, sizeof(facetT));
00960   if (qh facet_id == qh tracefacet_id)
00961     qh tracefacet= facet;
00962   facet->id= qh facet_id++;
00963   facet->neighbors= qh_setnew(qh hull_dim);
00964 #if !qh_COMPUTEfurthest
00965   facet->furthestdist= 0.0;
00966 #endif
00967 #if qh_MAXoutside
00968   if (qh FORCEoutput && qh APPROXhull)
00969     facet->maxoutside= qh MINoutside;
00970   else
00971     facet->maxoutside= qh DISTround;
00972 #endif
00973   facet->simplicial= True;
00974   facet->good= True;
00975   facet->newfacet= True;
00976   trace4((qh ferr, "qh_newfacet: created facet f%d\n", facet->id));
00977   return (facet);
00978 } /* newfacet */
00979 
00980 
00981 /*-<a                             href="qh-poly.htm#TOC"
00982   >-------------------------------</a><a name="newridge">-</a>
00983   
00984   qh_newridge()
00985     return a new ridge
00986 */
00987 ridgeT *qh_newridge(void) {
00988   ridgeT *ridge;
00989   void **freelistp;   /* used !qh_NOmem */
00990 
00991   qh_memalloc_(sizeof(ridgeT), freelistp, ridge, ridgeT);
00992   memset ((char *)ridge, 0, sizeof(ridgeT));
00993   zinc_(Ztotridges);
00994   if (qh ridge_id == 0xFFFFFF) {
00995     fprintf(qh ferr, "\
00996 qhull warning: more than %d ridges.  ID field overflows and two ridges\n\
00997 may have the same identifier.  Otherwise output ok.\n", 0xFFFFFF);
00998   }
00999   ridge->id= qh ridge_id++;     
01000   trace4((qh ferr, "qh_newridge: created ridge r%d\n", ridge->id));
01001   return (ridge);
01002 } /* newridge */
01003 
01004 
01005 /*-<a                             href="qh-poly.htm#TOC"
01006   >-------------------------------</a><a name="pointid">-</a>
01007   
01008   qh_pointid(  )
01009     return id for a point, 
01010     returns -3 if null, -2 if interior, or -1 if not known
01011 
01012   alternative code:
01013     unsigned long id;
01014     id= ((unsigned long)point - (unsigned long)qh.first_point)/qh.normal_size;
01015 
01016   notes:
01017     if point not in point array
01018       the code does a comparison of unrelated pointers.
01019 */
01020 int qh_pointid (pointT *point) {
01021   long offset, id;
01022 
01023   if (!point)
01024     id= -3;
01025   else if (point == qh interior_point)
01026     id= -2;
01027   else if (point >= qh first_point
01028   && point < qh first_point + qh num_points * qh hull_dim) {
01029     offset= point - qh first_point;
01030     id= offset / qh hull_dim;
01031   }else if ((id= qh_setindex (qh other_points, point)) != -1)
01032     id += qh num_points;
01033   else
01034     id= -1;
01035   return (int) id;
01036 } /* pointid */
01037   
01038 /*-<a                             href="qh-poly.htm#TOC"
01039   >-------------------------------</a><a name="removefacet">-</a>
01040   
01041   qh_removefacet( facet )
01042     unlinks facet from qh.facet_list,
01043 
01044   returns:
01045     updates qh.facet_list .newfacet_list .facet_next visible_list
01046     decrements qh.num_facets
01047 */
01048 void qh_removefacet(facetT *facet) {
01049   facetT *next= facet->next, *previous= facet->previous;
01050   
01051   if (facet == qh newfacet_list)
01052     qh newfacet_list= next;
01053   if (facet == qh facet_next)
01054     qh facet_next= next;
01055   if (facet == qh visible_list)
01056     qh visible_list= next; 
01057   if (previous) {
01058     previous->next= next;
01059     next->previous= previous;
01060   }else {  /* 1st facet in qh facet_list */
01061     qh facet_list= next;
01062     qh facet_list->previous= NULL;
01063   }
01064   qh num_facets--;
01065   trace4((qh ferr, "qh_removefacet: remove f%d from facet_list\n", facet->id));
01066 } /* removefacet */
01067 
01068 
01069 /*-<a                             href="qh-poly.htm#TOC"
01070   >-------------------------------</a><a name="removevertex">-</a>
01071   
01072   qh_removevertex( vertex )
01073     unlinks vertex from qh.vertex_list,
01074 
01075   returns:
01076     updates qh.vertex_list .newvertex_list 
01077     decrements qh.num_vertices
01078 */
01079 void qh_removevertex(vertexT *vertex) {
01080   vertexT *next= vertex->next, *previous= vertex->previous;
01081   
01082   if (vertex == qh newvertex_list)
01083     qh newvertex_list= next;
01084   if (previous) {
01085     previous->next= next;
01086     next->previous= previous;
01087   }else {  /* 1st vertex in qh vertex_list */
01088     qh vertex_list= vertex->next;
01089     qh vertex_list->previous= NULL;
01090   }
01091   qh num_vertices--;
01092   trace4((qh ferr, "qh_removevertex: remove v%d from vertex_list\n", vertex->id));
01093 } /* removevertex */
01094 
01095 
01096 /*-<a                             href="qh-poly.htm#TOC"
01097   >-------------------------------</a><a name="updatevertices">-</a>
01098   
01099   qh_updatevertices()
01100     update vertex neighbors and delete interior vertices
01101 
01102   returns:
01103     if qh.VERTEXneighbors
01104       updates neighbors for each vertex
01105     interior vertices added to qh.del_vertices for later partitioning
01106 
01107   design:
01108     if qh.VERTEXneighbors
01109       deletes references to visible facets from vertex neighbors
01110       appends new facets to the neighbor list for each vertex
01111       checks all vertices of visible facets
01112         removes visible facets from neighbor lists
01113         marks unused vertices for deletion
01114 */
01115 void qh_updatevertices (void) {
01116   facetT *newfacet= NULL, *neighbor, **neighborp, *visible;
01117   vertexT *vertex, **vertexp;
01118 
01119   trace3((qh ferr, "qh_updatevertices: delete interior vertices and update vertex->neighbors\n"));
01120   if (qh VERTEXneighbors) {
01121     FORALLvertex_(qh newvertex_list) {
01122       FOREACHneighbor_(vertex) {
01123         if (neighbor->visible) 
01124           SETref_(neighbor)= NULL;
01125       }
01126       qh_setcompact (vertex->neighbors);
01127     }
01128     FORALLnew_facets {
01129       FOREACHvertex_(newfacet->vertices)
01130         qh_setappend (&vertex->neighbors, newfacet);
01131     }
01132     FORALLvisible_facets {
01133       FOREACHvertex_(visible->vertices) {
01134         if (!vertex->newlist && !vertex->deleted) {
01135           FOREACHneighbor_(vertex) { /* this can happen under merging */
01136             if (!neighbor->visible)
01137               break;
01138           }
01139           if (neighbor)
01140             qh_setdel (vertex->neighbors, visible);
01141           else {
01142             vertex->deleted= True;
01143             qh_setappend (&qh del_vertices, vertex);
01144             trace2((qh ferr, "qh_updatevertices: delete vertex p%d (v%d) in f%d\n",
01145                   qh_pointid(vertex->point), vertex->id, visible->id));
01146           }
01147         }
01148       }
01149     }
01150   }else {  /* !VERTEXneighbors */
01151     FORALLvisible_facets {
01152       FOREACHvertex_(visible->vertices) {
01153         if (!vertex->newlist && !vertex->deleted) {
01154           vertex->deleted= True;
01155           qh_setappend (&qh del_vertices, vertex);
01156           trace2((qh ferr, "qh_updatevertices: delete vertex p%d (v%d) in f%d\n",
01157                   qh_pointid(vertex->point), vertex->id, visible->id));
01158         }
01159       }
01160     }
01161   }
01162 } /* updatevertices */
01163 
01164 
01165
 
Powered by Plone

This site conforms to the following standards: