Free cookie consent management tool by TermsFeed Policy Generator

source: branches/HeuristicLab.Analysis.AlgorithmBehavior/qhull-2012.1/src/libqhull/poly.c @ 10207

Last change on this file since 10207 was 10207, checked in by ascheibe, 11 years ago

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