ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/QuickHull/poly.c
Revision: 1138
Committed: Tue May 29 22:51:00 2007 UTC (17 years, 11 months ago) by chuckv
Content type: text/plain
File size: 37948 byte(s)
Log Message:
Addded QuickHull to cvs.

File Contents

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