Reverted incorrect merge (missing files)
[blender.git] / extern / qhull / src / poly.c
1 /*<html><pre>  -<a                             href="qh-poly.htm"
2   >-------------------------------</a><a name="TOP">-</a>
3
4    poly.c 
5    implements polygons and simplices
6
7    see qh-poly.htm, poly.h and qhull.h
8
9    infrequent code is in poly2.c 
10    (all but top 50 and their callers 12/3/95)
11
12    copyright (c) 1993-2002, The Geometry Center
13 */
14
15 #include "qhull_a.h"
16
17 /*======== functions in alphabetical order ==========*/
18
19 /*-<a                             href="qh-poly.htm#TOC"
20   >-------------------------------</a><a name="appendfacet">-</a>
21   
22   qh_appendfacet( facet )
23     appends facet to end of qh.facet_list,
24
25   returns:
26     updates qh.newfacet_list, facet_next, facet_list
27     increments qh.numfacets
28   
29   notes:
30     assumes qh.facet_list/facet_tail is defined (createsimplex)
31
32   see:
33     qh_removefacet()
34
35 */
36 void qh_appendfacet(facetT *facet) {
37   facetT *tail= qh facet_tail;
38
39   if (tail == qh newfacet_list)
40     qh newfacet_list= facet;
41   if (tail == qh facet_next)
42     qh facet_next= facet;
43   facet->previous= tail->previous;
44   facet->next= tail;
45   if (tail->previous)
46     tail->previous->next= facet;
47   else
48     qh facet_list= facet;
49   tail->previous= facet;
50   qh num_facets++;
51   trace4((qh ferr, "qh_appendfacet: append f%d to facet_list\n", facet->id));
52 } /* appendfacet */
53
54
55 /*-<a                             href="qh-poly.htm#TOC"
56   >-------------------------------</a><a name="appendvertex">-</a>
57   
58   qh_appendvertex( vertex )
59     appends vertex to end of qh.vertex_list,
60
61   returns:
62     sets vertex->newlist
63     updates qh.vertex_list, newvertex_list
64     increments qh.num_vertices
65
66   notes:
67     assumes qh.vertex_list/vertex_tail is defined (createsimplex)
68
69 */
70 void qh_appendvertex (vertexT *vertex) {
71   vertexT *tail= qh vertex_tail;
72
73   if (tail == qh newvertex_list)
74     qh newvertex_list= vertex;
75   vertex->newlist= True;
76   vertex->previous= tail->previous;
77   vertex->next= tail;
78   if (tail->previous)
79     tail->previous->next= vertex;
80   else
81     qh vertex_list= vertex;
82   tail->previous= vertex;
83   qh num_vertices++;
84   trace4((qh ferr, "qh_appendvertex: append v%d to vertex_list\n", vertex->id));
85 } /* appendvertex */
86
87
88 /*-<a                             href="qh-poly.htm#TOC"
89   >-------------------------------</a><a name="attachnewfacets">-</a>
90   
91   qh_attachnewfacets( )
92     attach horizon facets to new facets in qh.newfacet_list
93     newfacets have neighbor and ridge links to horizon but not vice versa
94     only needed for qh.ONLYgood
95
96   returns:
97     set qh.NEWfacets
98     horizon facets linked to new facets 
99       ridges changed from visible facets to new facets
100       simplicial ridges deleted
101     qh.visible_list, no ridges valid
102     facet->f.replace is a newfacet (if any)
103
104   design:
105     delete interior ridges and neighbor sets by
106       for each visible, non-simplicial facet
107         for each ridge
108           if last visit or if neighbor is simplicial
109             if horizon neighbor
110               delete ridge for horizon's ridge set
111             delete ridge
112         erase neighbor set
113     attach horizon facets and new facets by
114       for all new facets
115         if corresponding horizon facet is simplicial
116           locate corresponding visible facet {may be more than one}
117           link visible facet to new facet
118           replace visible facet with new facet in horizon
119         else it's non-simplicial
120           for all visible neighbors of the horizon facet
121             link visible neighbor to new facet
122             delete visible neighbor from horizon facet
123           append new facet to horizon's neighbors
124           the first ridge of the new facet is the horizon ridge
125           link the new facet into the horizon ridge
126 */
127 void qh_attachnewfacets (void ) {
128   facetT *newfacet= NULL, *neighbor, **neighborp, *horizon, *visible;
129   ridgeT *ridge, **ridgep;
130
131   qh NEWfacets= True;
132   trace3((qh ferr, "qh_attachnewfacets: delete interior ridges\n"));
133   qh visit_id++;
134   FORALLvisible_facets {
135     visible->visitid= qh visit_id;
136     if (visible->ridges) {
137       FOREACHridge_(visible->ridges) {
138         neighbor= otherfacet_(ridge, visible);
139         if (neighbor->visitid == qh visit_id
140             || (!neighbor->visible && neighbor->simplicial)) {
141           if (!neighbor->visible)  /* delete ridge for simplicial horizon */
142             qh_setdel (neighbor->ridges, ridge);
143           qh_setfree (&(ridge->vertices)); /* delete on 2nd visit */
144           qh_memfree (ridge, sizeof(ridgeT));
145         }
146       }
147       SETfirst_(visible->ridges)= NULL;
148     }
149     SETfirst_(visible->neighbors)= NULL;
150   }
151   trace1((qh ferr, "qh_attachnewfacets: attach horizon facets to new facets\n"));
152   FORALLnew_facets {
153     horizon= SETfirstt_(newfacet->neighbors, facetT);
154     if (horizon->simplicial) {
155       visible= NULL;
156       FOREACHneighbor_(horizon) {   /* may have more than one horizon ridge */
157         if (neighbor->visible) {
158           if (visible) {
159             if (qh_setequal_skip (newfacet->vertices, 0, horizon->vertices,
160                                   SETindex_(horizon->neighbors, neighbor))) {
161               visible= neighbor;
162               break;
163             }
164           }else
165             visible= neighbor;
166         }
167       }
168       if (visible) {
169         visible->f.replace= newfacet;
170         qh_setreplace (horizon->neighbors, visible, newfacet);
171       }else {
172         fprintf (qh ferr, "qhull internal error (qh_attachnewfacets): couldn't find visible facet for horizon f%d of newfacet f%d\n",
173                  horizon->id, newfacet->id);
174         qh_errexit2 (qh_ERRqhull, horizon, newfacet);
175       }
176     }else { /* non-simplicial, with a ridge for newfacet */
177       FOREACHneighbor_(horizon) {    /* may hold for many new facets */
178         if (neighbor->visible) {
179           neighbor->f.replace= newfacet;
180           qh_setdelnth (horizon->neighbors,
181                         SETindex_(horizon->neighbors, neighbor));
182           neighborp--; /* repeat */
183         }
184       }
185       qh_setappend (&horizon->neighbors, newfacet);
186       ridge= SETfirstt_(newfacet->ridges, ridgeT);
187       if (ridge->top == horizon)
188         ridge->bottom= newfacet;
189       else
190         ridge->top= newfacet;
191       }
192   } /* newfacets */
193   if (qh PRINTstatistics) {
194     FORALLvisible_facets {
195       if (!visible->f.replace) 
196         zinc_(Zinsidevisible);
197     }
198   }
199 } /* attachnewfacets */
200
201 /*-<a                             href="qh-poly.htm#TOC"
202   >-------------------------------</a><a name="checkflipped">-</a>
203   
204   qh_checkflipped( facet, dist, allerror )
205     checks facet orientation to interior point
206
207     if allerror set,
208       tests against qh.DISTround
209     else
210       tests against 0 since tested against DISTround before
211
212   returns:
213     False if it flipped orientation (sets facet->flipped)
214     distance if non-NULL
215 */
216 boolT qh_checkflipped (facetT *facet, realT *distp, boolT allerror) {
217   realT dist;
218
219   if (facet->flipped && !distp)
220     return False;
221   zzinc_(Zdistcheck);
222   qh_distplane(qh interior_point, facet, &dist);
223   if (distp)
224     *distp= dist;
225   if ((allerror && dist > -qh DISTround)|| (!allerror && dist >= 0.0)) {
226     facet->flipped= True;
227     zzinc_(Zflippedfacets);
228     trace0((qh ferr, "qh_checkflipped: facet f%d is flipped, distance= %6.12g during p%d\n",
229               facet->id, dist, qh furthest_id));
230     qh_precision ("flipped facet");
231     return False;
232   }
233   return True;
234 } /* checkflipped */
235
236 /*-<a                             href="qh-poly.htm#TOC"
237   >-------------------------------</a><a name="delfacet">-</a>
238   
239   qh_delfacet( facet )
240     removes facet from facet_list and frees up its memory
241
242   notes:
243     assumes vertices and ridges already freed
244 */
245 void qh_delfacet(facetT *facet) {
246   void **freelistp; /* used !qh_NOmem */
247
248   trace4((qh ferr, "qh_delfacet: delete f%d\n", facet->id));
249   if (facet == qh tracefacet)
250     qh tracefacet= NULL;
251   if (facet == qh GOODclosest)
252     qh GOODclosest= NULL;
253   qh_removefacet(facet);
254   if (!facet->tricoplanar || facet->keepcentrum) {
255     qh_memfree_(facet->normal, qh normal_size, freelistp);
256     if (qh CENTERtype == qh_ASvoronoi) {   /* uses macro calls */
257       qh_memfree_(facet->center, qh center_size, freelistp);
258     }else /* AScentrum */ {
259       qh_memfree_(facet->center, qh normal_size, freelistp);
260     }
261   }
262   qh_setfree(&(facet->neighbors));
263   if (facet->ridges)
264     qh_setfree(&(facet->ridges));
265   qh_setfree(&(facet->vertices));
266   if (facet->outsideset)
267     qh_setfree(&(facet->outsideset));
268   if (facet->coplanarset)
269     qh_setfree(&(facet->coplanarset));
270   qh_memfree_(facet, sizeof(facetT), freelistp);
271 } /* delfacet */
272
273
274 /*-<a                             href="qh-poly.htm#TOC"
275   >-------------------------------</a><a name="deletevisible">-</a>
276   
277   qh_deletevisible()
278     delete visible facets and vertices
279
280   returns:
281     deletes each facet and removes from facetlist
282     at exit, qh.visible_list empty (== qh.newfacet_list)
283
284   notes:
285     ridges already deleted
286     horizon facets do not reference facets on qh.visible_list
287     new facets in qh.newfacet_list
288     uses   qh.visit_id;
289 */
290 void qh_deletevisible (void /*qh visible_list*/) {
291   facetT *visible, *nextfacet;
292   vertexT *vertex, **vertexp;
293   int numvisible= 0, numdel= qh_setsize(qh del_vertices);
294
295   trace1((qh ferr, "qh_deletevisible: delete %d visible facets and %d vertices\n",
296          qh num_visible, numdel));
297   for (visible= qh visible_list; visible && visible->visible; 
298                 visible= nextfacet) { /* deleting current */
299     nextfacet= visible->next;        
300     numvisible++;
301     qh_delfacet(visible);
302   }
303   if (numvisible != qh num_visible) {
304     fprintf (qh ferr, "qhull internal error (qh_deletevisible): qh num_visible %d is not number of visible facets %d\n",
305              qh num_visible, numvisible);
306     qh_errexit (qh_ERRqhull, NULL, NULL);
307   }
308   qh num_visible= 0;
309   zadd_(Zvisfacettot, numvisible);
310   zmax_(Zvisfacetmax, numvisible);
311   zzadd_(Zdelvertextot, numdel);
312   zmax_(Zdelvertexmax, numdel);
313   FOREACHvertex_(qh del_vertices) 
314     qh_delvertex (vertex);
315   qh_settruncate (qh del_vertices, 0);
316 } /* deletevisible */
317
318 /*-<a                             href="qh-poly.htm#TOC"
319   >-------------------------------</a><a name="facetintersect">-</a>
320   
321   qh_facetintersect( facetA, facetB, skipa, skipB, prepend )
322     return vertices for intersection of two simplicial facets
323     may include 1 prepended entry (if more, need to settemppush)
324     
325   returns:
326     returns set of qh.hull_dim-1 + prepend vertices
327     returns skipped index for each test and checks for exactly one
328
329   notes:
330     does not need settemp since set in quick memory
331   
332   see also:
333     qh_vertexintersect and qh_vertexintersect_new
334     use qh_setnew_delnthsorted to get nth ridge (no skip information)
335
336   design:
337     locate skipped vertex by scanning facet A's neighbors
338     locate skipped vertex by scanning facet B's neighbors
339     intersect the vertex sets
340 */
341 setT *qh_facetintersect (facetT *facetA, facetT *facetB,
342                          int *skipA,int *skipB, int prepend) {
343   setT *intersect;
344   int dim= qh hull_dim, i, j;
345   facetT **neighborsA, **neighborsB;
346
347   neighborsA= SETaddr_(facetA->neighbors, facetT);
348   neighborsB= SETaddr_(facetB->neighbors, facetT);
349   i= j= 0;
350   if (facetB == *neighborsA++)
351     *skipA= 0;
352   else if (facetB == *neighborsA++)
353     *skipA= 1;
354   else if (facetB == *neighborsA++)
355     *skipA= 2;
356   else {
357     for (i= 3; i < dim; i++) {
358       if (facetB == *neighborsA++) {
359         *skipA= i;
360         break;
361       }
362     }
363   }
364   if (facetA == *neighborsB++)
365     *skipB= 0;
366   else if (facetA == *neighborsB++)
367     *skipB= 1;
368   else if (facetA == *neighborsB++)
369     *skipB= 2;
370   else {
371     for (j= 3; j < dim; j++) {
372       if (facetA == *neighborsB++) {
373         *skipB= j;
374         break;
375       }
376     }
377   }
378   if (i >= dim || j >= dim) {
379     fprintf (qh ferr, "qhull internal error (qh_facetintersect): f%d or f%d not in others neighbors\n",
380             facetA->id, facetB->id);
381     qh_errexit2 (qh_ERRqhull, facetA, facetB);
382   }
383   intersect= qh_setnew_delnthsorted (facetA->vertices, qh hull_dim, *skipA, prepend);
384   trace4((qh ferr, "qh_facetintersect: f%d skip %d matches f%d skip %d\n",
385           facetA->id, *skipA, facetB->id, *skipB));
386   return(intersect);
387 } /* facetintersect */
388
389 /*-<a                             href="qh-poly.htm#TOC"
390   >-------------------------------</a><a name="gethash">-</a>
391   
392   qh_gethash( hashsize, set, size, firstindex, skipelem )
393     return hashvalue for a set with firstindex and skipelem
394
395   notes:
396     assumes at least firstindex+1 elements
397     assumes skipelem is NULL, in set, or part of hash
398     
399     hashes memory addresses which may change over different runs of the same data
400     using sum for hash does badly in high d
401 */
402 unsigned qh_gethash (int hashsize, setT *set, int size, int firstindex, void *skipelem) {
403   void **elemp= SETelemaddr_(set, firstindex, void);
404   ptr_intT hash = 0, elem;
405   int i;
406
407   switch (size-firstindex) {
408   case 1:
409     hash= (ptr_intT)(*elemp) - (ptr_intT) skipelem;
410     break;
411   case 2:
412     hash= (ptr_intT)(*elemp) + (ptr_intT)elemp[1] - (ptr_intT) skipelem;
413     break;
414   case 3:
415     hash= (ptr_intT)(*elemp) + (ptr_intT)elemp[1] + (ptr_intT)elemp[2]
416       - (ptr_intT) skipelem;
417     break;
418   case 4:
419     hash= (ptr_intT)(*elemp) + (ptr_intT)elemp[1] + (ptr_intT)elemp[2]
420       + (ptr_intT)elemp[3] - (ptr_intT) skipelem;
421     break;
422   case 5:
423     hash= (ptr_intT)(*elemp) + (ptr_intT)elemp[1] + (ptr_intT)elemp[2]
424       + (ptr_intT)elemp[3] + (ptr_intT)elemp[4] - (ptr_intT) skipelem;
425     break;
426   case 6:
427     hash= (ptr_intT)(*elemp) + (ptr_intT)elemp[1] + (ptr_intT)elemp[2]
428       + (ptr_intT)elemp[3] + (ptr_intT)elemp[4]+ (ptr_intT)elemp[5]
429       - (ptr_intT) skipelem;
430     break;
431   default:
432     hash= 0;
433     i= 3;
434     do {     /* this is about 10% in 10-d */
435       if ((elem= (ptr_intT)*elemp++) != (ptr_intT)skipelem) {
436         hash ^= (elem << i) + (elem >> (32-i));
437         i += 3;
438         if (i >= 32)
439           i -= 32;
440       }
441     }while(*elemp);
442     break;
443   }
444   hash %= (ptr_intT) hashsize;
445   /* hash= 0; for debugging purposes */
446   return hash;
447 } /* gethash */
448
449 /*-<a                             href="qh-poly.htm#TOC"
450   >-------------------------------</a><a name="makenewfacet">-</a>
451   
452   qh_makenewfacet( vertices, toporient, horizon )
453     creates a toporient? facet from vertices
454
455   returns:
456     returns newfacet
457       adds newfacet to qh.facet_list
458       newfacet->vertices= vertices
459       if horizon
460         newfacet->neighbor= horizon, but not vice versa
461     newvertex_list updated with vertices
462 */
463 facetT *qh_makenewfacet(setT *vertices, boolT toporient,facetT *horizon) {
464   facetT *newfacet;
465   vertexT *vertex, **vertexp;
466
467   FOREACHvertex_(vertices) {
468     if (!vertex->newlist) {
469       qh_removevertex (vertex);
470       qh_appendvertex (vertex);
471     }
472   }
473   newfacet= qh_newfacet();
474   newfacet->vertices= vertices;
475   newfacet->toporient= toporient;
476   if (horizon)
477     qh_setappend(&(newfacet->neighbors), horizon);
478   qh_appendfacet(newfacet);
479   return(newfacet);
480 } /* makenewfacet */
481
482
483 /*-<a                             href="qh-poly.htm#TOC"
484   >-------------------------------</a><a name="makenewplanes">-</a>
485   
486   qh_makenewplanes()
487     make new hyperplanes for facets on qh.newfacet_list
488
489   returns:
490     all facets have hyperplanes or are marked for   merging
491     doesn't create hyperplane if horizon is coplanar (will merge)
492     updates qh.min_vertex if qh.JOGGLEmax
493
494   notes:
495     facet->f.samecycle is defined for facet->mergehorizon facets
496 */
497 void qh_makenewplanes (void /* newfacet_list */) {
498   facetT *newfacet;
499
500   FORALLnew_facets {
501     if (!newfacet->mergehorizon)
502       qh_setfacetplane (newfacet);  
503   }
504   if (qh JOGGLEmax < REALmax/2)  
505     minimize_(qh min_vertex, -wwval_(Wnewvertexmax));
506 } /* makenewplanes */
507
508 /*-<a                             href="qh-poly.htm#TOC"
509   >-------------------------------</a><a name="makenew_nonsimplicial">-</a>
510   
511   qh_makenew_nonsimplicial( visible, apex, numnew )
512     make new facets for ridges of a visible facet
513     
514   returns:
515     first newfacet, bumps numnew as needed
516     attaches new facets if !qh.ONLYgood
517     marks ridge neighbors for simplicial visible
518     if (qh.ONLYgood)
519       ridges on newfacet, horizon, and visible
520     else
521       ridge and neighbors between newfacet and   horizon
522       visible facet's ridges are deleted    
523
524   notes:
525     qh.visit_id if visible has already been processed
526     sets neighbor->seen for building f.samecycle
527       assumes all 'seen' flags initially false
528     
529   design:
530     for each ridge of visible facet
531       get neighbor of visible facet
532       if neighbor was already processed
533         delete the ridge (will delete all visible facets later)
534       if neighbor is a horizon facet
535         create a new facet
536         if neighbor coplanar
537           adds newfacet to f.samecycle for later merging
538         else 
539           updates neighbor's neighbor set
540           (checks for non-simplicial facet with multiple ridges to visible facet)
541         updates neighbor's ridge set
542         (checks for simplicial neighbor to non-simplicial visible facet)
543         (deletes ridge if neighbor is simplicial)
544           
545 */
546 #ifndef qh_NOmerge
547 facetT *qh_makenew_nonsimplicial (facetT *visible, vertexT *apex, int *numnew) {
548   void **freelistp; /* used !qh_NOmem */
549   ridgeT *ridge, **ridgep;
550   facetT *neighbor, *newfacet= NULL, *samecycle;
551   setT *vertices;
552   boolT toporient;
553   int ridgeid;
554
555   FOREACHridge_(visible->ridges) {
556     ridgeid= ridge->id;
557     neighbor= otherfacet_(ridge, visible);
558     if (neighbor->visible) {
559       if (!qh ONLYgood) {
560         if (neighbor->visitid == qh visit_id) {
561           qh_setfree (&(ridge->vertices));  /* delete on 2nd visit */
562           qh_memfree_(ridge, sizeof(ridgeT), freelistp);
563         }
564       }
565     }else {  /* neighbor is an horizon facet */
566       toporient= (ridge->top == visible);
567       vertices= qh_setnew (qh hull_dim); /* makes sure this is quick */
568       qh_setappend (&vertices, apex);
569       qh_setappend_set (&vertices, ridge->vertices);
570       newfacet= qh_makenewfacet(vertices, toporient, neighbor);
571       (*numnew)++;
572       if (neighbor->coplanar) {
573         newfacet->mergehorizon= True;
574         if (!neighbor->seen) {
575           newfacet->f.samecycle= newfacet;
576           neighbor->f.newcycle= newfacet;
577         }else {
578           samecycle= neighbor->f.newcycle;
579           newfacet->f.samecycle= samecycle->f.samecycle;
580           samecycle->f.samecycle= newfacet;
581         }
582       }
583       if (qh ONLYgood) {
584         if (!neighbor->simplicial)
585           qh_setappend(&(newfacet->ridges), ridge);
586       }else {  /* qh_attachnewfacets */
587         if (neighbor->seen) {
588           if (neighbor->simplicial) {
589             fprintf (qh ferr, "qhull internal error (qh_makenew_nonsimplicial): simplicial f%d sharing two ridges with f%d\n", 
590                    neighbor->id, visible->id);
591             qh_errexit2 (qh_ERRqhull, neighbor, visible);
592           }
593           qh_setappend (&(neighbor->neighbors), newfacet);
594         }else
595           qh_setreplace (neighbor->neighbors, visible, newfacet);
596         if (neighbor->simplicial) {
597           qh_setdel (neighbor->ridges, ridge);
598           qh_setfree (&(ridge->vertices)); 
599           qh_memfree (ridge, sizeof(ridgeT));
600         }else {
601           qh_setappend(&(newfacet->ridges), ridge);
602           if (toporient)
603             ridge->top= newfacet;
604           else
605             ridge->bottom= newfacet;
606         }
607       trace4((qh ferr, "qh_makenew_nonsimplicial: created facet f%d from v%d and r%d of horizon f%d\n",
608             newfacet->id, apex->id, ridgeid, neighbor->id));
609       }
610     }
611     neighbor->seen= True;        
612   } /* for each ridge */
613   if (!qh ONLYgood)
614     SETfirst_(visible->ridges)= NULL;
615   return newfacet;
616 } /* makenew_nonsimplicial */
617 #else /* qh_NOmerge */
618 facetT *qh_makenew_nonsimplicial (facetT *visible, vertexT *apex, int *numnew) {
619   return NULL;
620 }
621 #endif /* qh_NOmerge */
622
623 /*-<a                             href="qh-poly.htm#TOC"
624   >-------------------------------</a><a name="makenew_simplicial">-</a>
625   
626   qh_makenew_simplicial( visible, apex, numnew )
627     make new facets for simplicial visible facet and apex
628
629   returns:
630     attaches new facets if (!qh.ONLYgood)
631       neighbors between newfacet and horizon
632
633   notes:
634     nop if neighbor->seen or neighbor->visible (see qh_makenew_nonsimplicial)
635
636   design:
637     locate neighboring horizon facet for visible facet
638     determine vertices and orientation
639     create new facet
640     if coplanar,
641       add new facet to f.samecycle
642     update horizon facet's neighbor list        
643 */
644 facetT *qh_makenew_simplicial (facetT *visible, vertexT *apex, int *numnew) {
645   facetT *neighbor, **neighborp, *newfacet= NULL;
646   setT *vertices;
647   boolT flip, toporient;
648   int horizonskip, visibleskip;
649
650   FOREACHneighbor_(visible) {
651     if (!neighbor->seen && !neighbor->visible) {
652       vertices= qh_facetintersect(neighbor,visible, &horizonskip, &visibleskip, 1);
653       SETfirst_(vertices)= apex;
654       flip= ((horizonskip & 0x1) ^ (visibleskip & 0x1));
655       if (neighbor->toporient)         
656         toporient= horizonskip & 0x1;
657       else
658         toporient= (horizonskip & 0x1) ^ 0x1;
659       newfacet= qh_makenewfacet(vertices, toporient, neighbor);
660       (*numnew)++;
661       if (neighbor->coplanar && (qh PREmerge || qh MERGEexact)) {
662 #ifndef qh_NOmerge
663         newfacet->f.samecycle= newfacet;
664         newfacet->mergehorizon= True;
665 #endif
666       }
667       if (!qh ONLYgood)
668         SETelem_(neighbor->neighbors, horizonskip)= newfacet;
669       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",
670             newfacet->id, toporient, apex->id, neighbor->id, horizonskip,
671               neighbor->toporient, visible->id, visibleskip, flip));
672     }
673   }
674   return newfacet;
675 } /* makenew_simplicial */
676
677 /*-<a                             href="qh-poly.htm#TOC"
678   >-------------------------------</a><a name="matchneighbor">-</a>
679   
680   qh_matchneighbor( newfacet, newskip, hashsize, hashcount )
681     either match subridge of newfacet with neighbor or add to hash_table
682
683   returns:
684     duplicate ridges are unmatched and marked by qh_DUPLICATEridge
685
686   notes:
687     ridge is newfacet->vertices w/o newskip vertex
688     do not allocate memory (need to free hash_table cleanly)
689     uses linear hash chains
690   
691   see also:
692     qh_matchduplicates
693
694   design:
695     for each possible matching facet in qh.hash_table
696       if vertices match
697         set ismatch, if facets have opposite orientation
698         if ismatch and matching facet doesn't have a match
699           match the facets by updating their neighbor sets
700         else
701           indicate a duplicate ridge
702           set facet hyperplane for later testing
703           add facet to hashtable
704           unless the other facet was already a duplicate ridge
705             mark both facets with a duplicate ridge
706             add other facet (if defined) to hash table
707 */
708 void qh_matchneighbor (facetT *newfacet, int newskip, int hashsize, int *hashcount) {
709   boolT newfound= False;   /* True, if new facet is already in hash chain */
710   boolT same, ismatch;
711   int hash, scan;
712   facetT *facet, *matchfacet;
713   int skip, matchskip;
714
715   hash= (int)qh_gethash (hashsize, newfacet->vertices, qh hull_dim, 1, 
716                      SETelem_(newfacet->vertices, newskip));
717   trace4((qh ferr, "qh_matchneighbor: newfacet f%d skip %d hash %d hashcount %d\n",
718           newfacet->id, newskip, hash, *hashcount));
719   zinc_(Zhashlookup);
720   for (scan= hash; (facet= SETelemt_(qh hash_table, scan, facetT)); 
721        scan= (++scan >= hashsize ? 0 : scan)) {
722     if (facet == newfacet) {
723       newfound= True;
724       continue;
725     }
726     zinc_(Zhashtests);
727     if (qh_matchvertices (1, newfacet->vertices, newskip, facet->vertices, &skip, &same)) {
728       if (SETelem_(newfacet->vertices, newskip) == 
729           SETelem_(facet->vertices, skip)) {
730         qh_precision ("two facets with the same vertices");
731         fprintf (qh ferr, "qhull precision error: Vertex sets are the same for f%d and f%d.  Can not force output.\n",
732           facet->id, newfacet->id);
733         qh_errexit2 (qh_ERRprec, facet, newfacet);
734       }
735       ismatch= (same == (newfacet->toporient ^ facet->toporient));
736       matchfacet= SETelemt_(facet->neighbors, skip, facetT);
737       if (ismatch && !matchfacet) {
738         SETelem_(facet->neighbors, skip)= newfacet;
739         SETelem_(newfacet->neighbors, newskip)= facet;
740         (*hashcount)--;
741         trace4((qh ferr, "qh_matchneighbor: f%d skip %d matched with new f%d skip %d\n",
742            facet->id, skip, newfacet->id, newskip));
743         return;
744       }
745       if (!qh PREmerge && !qh MERGEexact) {
746         qh_precision ("a ridge with more than two neighbors");
747         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",
748                  facet->id, newfacet->id, getid_(matchfacet));
749         qh_errexit2 (qh_ERRprec, facet, newfacet);
750       }
751       SETelem_(newfacet->neighbors, newskip)= qh_DUPLICATEridge;
752       newfacet->dupridge= True;
753       if (!newfacet->normal)
754         qh_setfacetplane (newfacet);
755       qh_addhash (newfacet, qh hash_table, hashsize, hash);
756       (*hashcount)++;
757       if (!facet->normal)
758         qh_setfacetplane (facet);
759       if (matchfacet != qh_DUPLICATEridge) {
760         SETelem_(facet->neighbors, skip)= qh_DUPLICATEridge;
761         facet->dupridge= True;
762         if (!facet->normal)
763           qh_setfacetplane (facet);
764         if (matchfacet) {
765           matchskip= qh_setindex (matchfacet->neighbors, facet);
766           SETelem_(matchfacet->neighbors, matchskip)= qh_DUPLICATEridge;
767           matchfacet->dupridge= True;
768           if (!matchfacet->normal)
769             qh_setfacetplane (matchfacet);
770           qh_addhash (matchfacet, qh hash_table, hashsize, hash);
771           *hashcount += 2;
772         }
773       }
774       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",
775            newfacet->id, newskip, facet->id, skip, 
776            (matchfacet == qh_DUPLICATEridge ? -2 : getid_(matchfacet)), 
777            ismatch, hash));
778       return; /* end of duplicate ridge */
779     }
780   }
781   if (!newfound) 
782     SETelem_(qh hash_table, scan)= newfacet;  /* same as qh_addhash */
783   (*hashcount)++;
784   trace4((qh ferr, "qh_matchneighbor: no match for f%d skip %d at hash %d\n",
785            newfacet->id, newskip, hash));
786 } /* matchneighbor */
787
788
789 /*-<a                             href="qh-poly.htm#TOC"
790   >-------------------------------</a><a name="matchnewfacets">-</a>
791   
792   qh_matchnewfacets()
793     match newfacets in qh.newfacet_list to their newfacet neighbors
794
795   returns:
796     qh.newfacet_list with full neighbor sets
797       get vertices with nth neighbor by deleting nth vertex
798     if qh.PREmerge/MERGEexact or qh.FORCEoutput 
799       sets facet->flippped if flipped normal (also prevents point partitioning)
800     if duplicate ridges and qh.PREmerge/MERGEexact
801       sets facet->dupridge
802       missing neighbor links identifies extra ridges to be merging (qh_MERGEridge)
803
804   notes:
805     newfacets already have neighbor[0] (horizon facet)
806     assumes qh.hash_table is NULL
807     vertex->neighbors has not been updated yet
808     do not allocate memory after qh.hash_table (need to free it cleanly)
809
810   design:
811     delete neighbor sets for all new facets
812     initialize a hash table
813     for all new facets
814       match facet with neighbors
815     if unmatched facets (due to duplicate ridges)
816       for each new facet with a duplicate ridge
817         match it with a facet
818     check for flipped facets
819 */
820 void qh_matchnewfacets (void /* qh newfacet_list */) {
821   int numnew=0, hashcount=0, newskip;
822   facetT *newfacet, *neighbor;
823   int dim= qh hull_dim, hashsize, neighbor_i, neighbor_n;
824   setT *neighbors;
825 #ifndef qh_NOtrace
826   int facet_i, facet_n, numfree= 0;
827   facetT *facet;
828 #endif
829   
830   trace1((qh ferr, "qh_matchnewfacets: match neighbors for new facets.\n"));
831   FORALLnew_facets {
832     numnew++;
833     {  /* inline qh_setzero (newfacet->neighbors, 1, qh hull_dim); */
834       neighbors= newfacet->neighbors;
835       neighbors->e[neighbors->maxsize].i= dim+1; /*may be overwritten*/
836       memset ((char *)SETelemaddr_(neighbors, 1, void), 0, dim * SETelemsize);
837     }    
838   }
839   qh_newhashtable (numnew*(qh hull_dim-1)); /* twice what is normally needed,
840                                      but every ridge could be DUPLICATEridge */
841   hashsize= qh_setsize (qh hash_table);
842   FORALLnew_facets {
843     for (newskip=1; newskip<qh hull_dim; newskip++) /* furthest/horizon already matched */
844       qh_matchneighbor (newfacet, newskip, hashsize, &hashcount);
845 #if 0   /* use the following to trap hashcount errors */
846     {
847       int count= 0, k;
848       facetT *facet, *neighbor;
849
850       count= 0;
851       FORALLfacet_(qh newfacet_list) {  /* newfacet already in use */
852         for (k=1; k < qh hull_dim; k++) {
853           neighbor= SETelemt_(facet->neighbors, k, facetT);
854           if (!neighbor || neighbor == qh_DUPLICATEridge)
855             count++;
856         }
857         if (facet == newfacet)
858           break;
859       }
860       if (count != hashcount) {
861         fprintf (qh ferr, "qh_matchnewfacets: after adding facet %d, hashcount %d != count %d\n",
862                  newfacet->id, hashcount, count);
863         qh_errexit (qh_ERRqhull, newfacet, NULL);
864       }
865     }
866 #endif  /* end of trap code */
867   }
868   if (hashcount) {
869     FORALLnew_facets {
870       if (newfacet->dupridge) {
871         FOREACHneighbor_i_(newfacet) {
872           if (neighbor == qh_DUPLICATEridge) {
873             qh_matchduplicates (newfacet, neighbor_i, hashsize, &hashcount);
874                     /* this may report MERGEfacet */
875           }
876         }
877       }
878     }
879   }
880   if (hashcount) {
881     fprintf (qh ferr, "qhull internal error (qh_matchnewfacets): %d neighbors did not match up\n",
882         hashcount);
883     qh_printhashtable (qh ferr);
884     qh_errexit (qh_ERRqhull, NULL, NULL);
885   }
886 #ifndef qh_NOtrace
887   if (qh IStracing >= 2) {
888     FOREACHfacet_i_(qh hash_table) {
889       if (!facet)
890         numfree++;
891     }
892     fprintf (qh ferr, "qh_matchnewfacets: %d new facets, %d unused hash entries .  hashsize %d\n",
893              numnew, numfree, qh_setsize (qh hash_table));
894   }
895 #endif /* !qh_NOtrace */
896   qh_setfree (&qh hash_table);
897   if (qh PREmerge || qh MERGEexact) {
898     if (qh IStracing >= 4)
899       qh_printfacetlist (qh newfacet_list, NULL, qh_ALL);
900     FORALLnew_facets {
901       if (newfacet->normal)
902         qh_checkflipped (newfacet, NULL, qh_ALL);
903     }
904   }else if (qh FORCEoutput)
905     qh_checkflipped_all (qh newfacet_list);  /* prints warnings for flipped */
906 } /* matchnewfacets */
907
908     
909 /*-<a                             href="qh-poly.htm#TOC"
910   >-------------------------------</a><a name="matchvertices">-</a>
911   
912   qh_matchvertices( firstindex, verticesA, skipA, verticesB, skipB, same )
913     tests whether vertices match with a single skip
914     starts match at firstindex since all new facets have a common vertex
915
916   returns:
917     true if matched vertices
918     skip index for each set
919     sets same iff vertices have the same orientation
920
921   notes:
922     assumes skipA is in A and both sets are the same size
923
924   design:
925     set up pointers
926     scan both sets checking for a match
927     test orientation
928 */
929 boolT qh_matchvertices (int firstindex, setT *verticesA, int skipA, 
930        setT *verticesB, int *skipB, boolT *same) {
931   vertexT **elemAp, **elemBp, **skipBp=NULL, **skipAp;
932
933   elemAp= SETelemaddr_(verticesA, firstindex, vertexT);
934   elemBp= SETelemaddr_(verticesB, firstindex, vertexT);
935   skipAp= SETelemaddr_(verticesA, skipA, vertexT);
936   do if (elemAp != skipAp) {
937     while (*elemAp != *elemBp++) {
938       if (skipBp)
939         return False;
940       skipBp= elemBp;  /* one extra like FOREACH */
941     }
942   }while(*(++elemAp));
943   if (!skipBp)
944     skipBp= ++elemBp;
945   *skipB= SETindex_(verticesB, skipB);
946   *same= !(((ptr_intT)skipA & 0x1) ^ ((ptr_intT)*skipB & 0x1));
947   trace4((qh ferr, "qh_matchvertices: matched by skip %d (v%d) and skip %d (v%d) same? %d\n",
948           skipA, (*skipAp)->id, *skipB, (*(skipBp-1))->id, *same));
949   return (True);
950 } /* matchvertices */
951
952 /*-<a                             href="qh-poly.htm#TOC"
953   >-------------------------------</a><a name="newfacet">-</a>
954   
955   qh_newfacet()
956     return a new facet 
957
958   returns:
959     all fields initialized or cleared   (NULL)
960     preallocates neighbors set
961 */
962 facetT *qh_newfacet(void) {
963   facetT *facet;
964   void **freelistp; /* used !qh_NOmem */
965   
966   qh_memalloc_(sizeof(facetT), freelistp, facet, facetT);
967   memset ((char *)facet, 0, sizeof(facetT));
968   if (qh facet_id == qh tracefacet_id)
969     qh tracefacet= facet;
970   facet->id= qh facet_id++;
971   facet->neighbors= qh_setnew(qh hull_dim);
972 #if !qh_COMPUTEfurthest
973   facet->furthestdist= 0.0;
974 #endif
975 #if qh_MAXoutside
976   if (qh FORCEoutput && qh APPROXhull)
977     facet->maxoutside= qh MINoutside;
978   else
979     facet->maxoutside= qh DISTround;
980 #endif
981   facet->simplicial= True;
982   facet->good= True;
983   facet->newfacet= True;
984   trace4((qh ferr, "qh_newfacet: created facet f%d\n", facet->id));
985   return (facet);
986 } /* newfacet */
987
988
989 /*-<a                             href="qh-poly.htm#TOC"
990   >-------------------------------</a><a name="newridge">-</a>
991   
992   qh_newridge()
993     return a new ridge
994 */
995 ridgeT *qh_newridge(void) {
996   ridgeT *ridge;
997   void **freelistp;   /* used !qh_NOmem */
998
999   qh_memalloc_(sizeof(ridgeT), freelistp, ridge, ridgeT);
1000   memset ((char *)ridge, 0, sizeof(ridgeT));
1001   zinc_(Ztotridges);
1002   if (qh ridge_id == 0xFFFFFF) {
1003     fprintf(qh ferr, "\
1004 qhull warning: more than %d ridges.  ID field overflows and two ridges\n\
1005 may have the same identifier.  Otherwise output ok.\n", 0xFFFFFF);
1006   }
1007   ridge->id= qh ridge_id++;     
1008   trace4((qh ferr, "qh_newridge: created ridge r%d\n", ridge->id));
1009   return (ridge);
1010 } /* newridge */
1011
1012
1013 /*-<a                             href="qh-poly.htm#TOC"
1014   >-------------------------------</a><a name="pointid">-</a>
1015   
1016   qh_pointid(  )
1017     return id for a point, 
1018     returns -3 if null, -2 if interior, or -1 if not known
1019
1020   alternative code:
1021     unsigned long id;
1022     id= ((unsigned long)point - (unsigned long)qh.first_point)/qh.normal_size;
1023
1024   notes:
1025     if point not in point array
1026       the code does a comparison of unrelated pointers.
1027 */
1028 int qh_pointid (pointT *point) {
1029   long offset, id;
1030
1031   if (!point)
1032     id= -3;
1033   else if (point == qh interior_point)
1034     id= -2;
1035   else if (point >= qh first_point
1036   && point < qh first_point + qh num_points * qh hull_dim) {
1037     offset= point - qh first_point;
1038     id= offset / qh hull_dim;
1039   }else if ((id= qh_setindex (qh other_points, point)) != -1)
1040     id += qh num_points;
1041   else
1042     id= -1;
1043   return (int) id;
1044 } /* pointid */
1045   
1046 /*-<a                             href="qh-poly.htm#TOC"
1047   >-------------------------------</a><a name="removefacet">-</a>
1048   
1049   qh_removefacet( facet )
1050     unlinks facet from qh.facet_list,
1051
1052   returns:
1053     updates qh.facet_list .newfacet_list .facet_next visible_list
1054     decrements qh.num_facets
1055
1056   see:
1057     qh_appendfacet
1058 */
1059 void qh_removefacet(facetT *facet) {
1060   facetT *next= facet->next, *previous= facet->previous;
1061   
1062   if (facet == qh newfacet_list)
1063     qh newfacet_list= next;
1064   if (facet == qh facet_next)
1065     qh facet_next= next;
1066   if (facet == qh visible_list)
1067     qh visible_list= next; 
1068   if (previous) {
1069     previous->next= next;
1070     next->previous= previous;
1071   }else {  /* 1st facet in qh facet_list */
1072     qh facet_list= next;
1073     qh facet_list->previous= NULL;
1074   }
1075   qh num_facets--;
1076   trace4((qh ferr, "qh_removefacet: remove f%d from facet_list\n", facet->id));
1077 } /* removefacet */
1078
1079
1080 /*-<a                             href="qh-poly.htm#TOC"
1081   >-------------------------------</a><a name="removevertex">-</a>
1082   
1083   qh_removevertex( vertex )
1084     unlinks vertex from qh.vertex_list,
1085
1086   returns:
1087     updates qh.vertex_list .newvertex_list 
1088     decrements qh.num_vertices
1089 */
1090 void qh_removevertex(vertexT *vertex) {
1091   vertexT *next= vertex->next, *previous= vertex->previous;
1092   
1093   if (vertex == qh newvertex_list)
1094     qh newvertex_list= next;
1095   if (previous) {
1096     previous->next= next;
1097     next->previous= previous;
1098   }else {  /* 1st vertex in qh vertex_list */
1099     qh vertex_list= vertex->next;
1100     qh vertex_list->previous= NULL;
1101   }
1102   qh num_vertices--;
1103   trace4((qh ferr, "qh_removevertex: remove v%d from vertex_list\n", vertex->id));
1104 } /* removevertex */
1105
1106
1107 /*-<a                             href="qh-poly.htm#TOC"
1108   >-------------------------------</a><a name="updatevertices">-</a>
1109   
1110   qh_updatevertices()
1111     update vertex neighbors and delete interior vertices
1112
1113   returns:
1114     if qh.VERTEXneighbors, updates neighbors for each vertex
1115       if qh.newvertex_list, 
1116          removes visible neighbors  from vertex neighbors
1117       if qh.newfacet_list
1118          adds new facets to vertex neighbors
1119     if qh.visible_list
1120        interior vertices added to qh.del_vertices for later partitioning
1121
1122   design:
1123     if qh.VERTEXneighbors
1124       deletes references to visible facets from vertex neighbors
1125       appends new facets to the neighbor list for each vertex
1126       checks all vertices of visible facets
1127         removes visible facets from neighbor lists
1128         marks unused vertices for deletion
1129 */
1130 void qh_updatevertices (void /*qh newvertex_list, newfacet_list, visible_list*/) {
1131   facetT *newfacet= NULL, *neighbor, **neighborp, *visible;
1132   vertexT *vertex, **vertexp;
1133
1134   trace3((qh ferr, "qh_updatevertices: delete interior vertices and update vertex->neighbors\n"));
1135   if (qh VERTEXneighbors) {
1136     FORALLvertex_(qh newvertex_list) {
1137       FOREACHneighbor_(vertex) {
1138         if (neighbor->visible) 
1139           SETref_(neighbor)= NULL;
1140       }
1141       qh_setcompact (vertex->neighbors);
1142     }
1143     FORALLnew_facets {
1144       FOREACHvertex_(newfacet->vertices)
1145         qh_setappend (&vertex->neighbors, newfacet);
1146     }
1147     FORALLvisible_facets {
1148       FOREACHvertex_(visible->vertices) {
1149         if (!vertex->newlist && !vertex->deleted) {
1150           FOREACHneighbor_(vertex) { /* this can happen under merging */
1151             if (!neighbor->visible)
1152               break;
1153           }
1154           if (neighbor)
1155             qh_setdel (vertex->neighbors, visible);
1156           else {
1157             vertex->deleted= True;
1158             qh_setappend (&qh del_vertices, vertex);
1159             trace2((qh ferr, "qh_updatevertices: delete vertex p%d (v%d) in f%d\n",
1160                   qh_pointid(vertex->point), vertex->id, visible->id));
1161           }
1162         }
1163       }
1164     }
1165   }else {  /* !VERTEXneighbors */
1166     FORALLvisible_facets {
1167       FOREACHvertex_(visible->vertices) {
1168         if (!vertex->newlist && !vertex->deleted) {
1169           vertex->deleted= True;
1170           qh_setappend (&qh del_vertices, vertex);
1171           trace2((qh ferr, "qh_updatevertices: delete vertex p%d (v%d) in f%d\n",
1172                   qh_pointid(vertex->point), vertex->id, visible->id));
1173         }
1174       }
1175     }
1176   }
1177 } /* updatevertices */
1178
1179
1180