/*
--------------------------------- merge.c merges non-convex facets see qh-merge.htm and merge.h other modules call qh_premerge() and qh_postmerge() the user may call qh_postmerge() to perform additional merges. To remove deleted facets and vertices (qhull() in libqhull.c): qh_partitionvisible(!qh_ALL, &numoutside); // visible_list, newfacet_list qh_deletevisible(); // qh.visible_list qh_resetlists(False, qh_RESETvisible); // qh.visible_list newvertex_list newfacet_list assumes qh.CENTERtype= centrum merges occur in qh_mergefacet and in qh_mergecycle vertex->neighbors not set until the first merge occurs Copyright (c) 1993-2020 C.B. Barber. $Id: //main/2019/qhull/src/libqhull/merge.c#12 $$Change: 2958 $ $DateTime: 2020/05/26 16:17:49 $$Author: bbarber $ */ #include "qhull_a.h" #ifndef qh_NOmerge /* MRGnone, etc. */ const char *mergetypes[]= { "none", "coplanar", "anglecoplanar", "concave", "concavecoplanar", "twisted", "flip", "dupridge", "subridge", "vertices", "degen", "redundant", "mirror", "coplanarhorizon", }; /*===== functions(alphabetical after premerge and postmerge) ======*/ /*--------------------------------- qh_premerge( apexpointid, maxcentrum ) pre-merge nonconvex facets in qh.newfacet_list for apexpointid maxcentrum defines coplanar and concave (qh_test_appendmerge) returns: deleted facets added to qh.visible_list with facet->visible set notes: only called by qh_addpoint uses globals, qh.MERGEexact, qh.PREmerge design: mark dupridges in qh.newfacet_list merge facet cycles in qh.newfacet_list merge dupridges and concave facets in qh.newfacet_list check merged facet cycles for degenerate and redundant facets merge degenerate and redundant facets collect coplanar and concave facets merge concave, coplanar, degenerate, and redundant facets */ void qh_premerge(int apexpointid, realT maxcentrum, realT maxangle /* qh.newfacet_list */) { boolT othermerge= False; if (qh ZEROcentrum && qh_checkzero(!qh_ALL)) return; trace2((qh ferr, 2008, "qh_premerge: premerge centrum %2.2g angle %4.4g for apex p%d newfacet_list f%d\n", maxcentrum, maxangle, apexpointid, getid_(qh newfacet_list))); if (qh IStracing >= 4 && qh num_facets < 100) qh_printlists(); qh centrum_radius= maxcentrum; qh cos_max= maxangle; if (qh hull_dim >=3) { qh_mark_dupridges(qh newfacet_list, qh_ALL); /* facet_mergeset */ qh_mergecycle_all(qh newfacet_list, &othermerge); qh_forcedmerges(&othermerge /* qh.facet_mergeset */); }else /* qh.hull_dim == 2 */ qh_mergecycle_all(qh newfacet_list, &othermerge); qh_flippedmerges(qh newfacet_list, &othermerge); if (!qh MERGEexact || zzval_(Ztotmerge)) { zinc_(Zpremergetot); qh POSTmerging= False; qh_getmergeset_initial(qh newfacet_list); qh_all_merges(othermerge, False); } } /* premerge */ /*--------------------------------- qh_postmerge( reason, maxcentrum, maxangle, vneighbors ) post-merge nonconvex facets as defined by maxcentrum and maxangle 'reason' is for reporting progress if vneighbors ('Qv'), calls qh_test_vneighbors at end of qh_all_merge from qh_postmerge returns: if first call (qh.visible_list != qh.facet_list), builds qh.facet_newlist, qh.newvertex_list deleted facets added to qh.visible_list with facet->visible qh.visible_list == qh.facet_list notes: called by qh_qhull after qh_buildhull called if a merge may be needed due to qh.MERGEexact ('Qx'), qh_DIMreduceBuild, POSTmerge (e.g., 'Cn'), or TESTvneighbors ('Qv') if firstmerge, calls qh_reducevertices before qh_getmergeset design: if first call set qh.visible_list and qh.newfacet_list to qh.facet_list add all facets to qh.newfacet_list mark non-simplicial facets, facet->newmerge set qh.newvertext_list to qh.vertex_list add all vertices to qh.newvertex_list if a pre-merge occurred set vertex->delridge {will retest the ridge} if qh.MERGEexact call qh_reducevertices() if no pre-merging merge flipped facets determine non-convex facets merge all non-convex facets */ void qh_postmerge(const char *reason, realT maxcentrum, realT maxangle, boolT vneighbors) { facetT *newfacet; boolT othermerges= False; vertexT *vertex; if (qh REPORTfreq || qh IStracing) { qh_buildtracing(NULL, NULL); qh_printsummary(qh ferr); if (qh PRINTstatistics) qh_printallstatistics(qh ferr, "reason"); qh_fprintf(qh ferr, 8062, "\n%s with 'C%.2g' and 'A%.2g'\n", reason, maxcentrum, maxangle); } trace2((qh ferr, 2009, "qh_postmerge: postmerge. test vneighbors? %d\n", vneighbors)); qh centrum_radius= maxcentrum; qh cos_max= maxangle; qh POSTmerging= True; if (qh visible_list != qh facet_list) { /* first call due to qh_buildhull, multiple calls if qh.POSTmerge */ qh NEWfacets= True; qh visible_list= qh newfacet_list= qh facet_list; FORALLnew_facets { /* all facets are new facets for qh_postmerge */ newfacet->newfacet= True; if (!newfacet->simplicial) newfacet->newmerge= True; /* test f.vertices for 'delridge'. 'newmerge' was cleared at end of qh_all_merges */ zinc_(Zpostfacets); } qh newvertex_list= qh vertex_list; FORALLvertices vertex->newfacet= True; if (qh VERTEXneighbors) { /* a merge has occurred */ if (qh MERGEexact && qh hull_dim <= qh_DIMreduceBuild) qh_reducevertices(); /* qh_all_merges did not call qh_reducevertices for v.delridge */ } if (!qh PREmerge && !qh MERGEexact) qh_flippedmerges(qh newfacet_list, &othermerges); } qh_getmergeset_initial(qh newfacet_list); qh_all_merges(False, vneighbors); /* calls qh_reducevertices before exiting */ FORALLnew_facets newfacet->newmerge= False; /* Was True if no vertex in f.vertices was 'delridge' */ } /* post_merge */ /*--------------------------------- qh_all_merges( othermerge, vneighbors ) merge all non-convex facets set othermerge if already merged facets (calls qh_reducevertices) if vneighbors ('Qv' at qh.POSTmerge) tests vertex neighbors for convexity at end (qh_test_vneighbors) qh.facet_mergeset lists the non-convex ridges in qh_newfacet_list qh.degen_mergeset is defined if qh.MERGEexact && !qh.POSTmerging, does not merge coplanar facets returns: deleted facets added to qh.visible_list with facet->visible deleted vertices added qh.delvertex_list with vertex->delvertex notes: unless !qh.MERGEindependent, merges facets in independent sets uses qh.newfacet_list as implicit argument since merges call qh_removefacet() [apr'19] restored qh_setdellast in place of qh_next_facetmerge. Much faster for post-merge design: while merges occur for each merge in qh.facet_mergeset unless one of the facets was already merged in this pass merge the facets test merged facets for additional merges add merges to qh.facet_mergeset if qh.POSTmerging periodically call qh_reducevertices to reduce extra vertices and redundant vertices after each pass, if qh.VERTEXneighbors if qh.POSTmerging or was a merge with qh.hull_dim<=5 call qh_reducevertices update qh.facet_mergeset if degenredundant merges if 'Qv' and qh.POSTmerging test vertex neighbors for convexity */ void qh_all_merges(boolT othermerge, boolT vneighbors) { facetT *facet1, *facet2, *newfacet; mergeT *merge; boolT wasmerge= False, isreduce; void **freelistp; /* used if !qh_NOmem by qh_memfree_() */ vertexT *vertex; realT angle, distance; mergeType mergetype; int numcoplanar=0, numconcave=0, numconcavecoplanar= 0, numdegenredun= 0, numnewmerges= 0, numtwisted= 0; trace2((qh ferr, 2010, "qh_all_merges: starting to merge %d facet and %d degenerate merges for new facets f%d, othermerge? %d\n", qh_setsize(qh facet_mergeset), qh_setsize(qh degen_mergeset), getid_(qh newfacet_list), othermerge)); while (True) { wasmerge= False; while (qh_setsize(qh facet_mergeset) > 0 || qh_setsize(qh degen_mergeset) > 0) { if (qh_setsize(qh degen_mergeset) > 0) { numdegenredun += qh_merge_degenredundant(); wasmerge= True; } while ((merge= (mergeT *)qh_setdellast(qh facet_mergeset))) { facet1= merge->facet1; facet2= merge->facet2; vertex= merge->vertex1; /* not used for qh.facet_mergeset*/ mergetype= merge->mergetype; angle= merge->angle; distance= merge->distance; qh_memfree_(merge, (int)sizeof(mergeT), freelistp); /* 'merge' is invalid */ if (facet1->visible || facet2->visible) { trace3((qh ferr, 3045, "qh_all_merges: drop merge of f%d (del? %d) into f%d (del? %d) mergetype %d, dist %4.4g, angle %4.4g. One or both facets is deleted\n", facet1->id, facet1->visible, facet2->id, facet2->visible, mergetype, distance, angle)); continue; }else if (mergetype == MRGcoplanar || mergetype == MRGanglecoplanar) { if (qh MERGEindependent) { if ((!facet1->tested && facet1->newfacet) || (!facet2->tested && facet2->newfacet)) { trace3((qh ferr, 3064, "qh_all_merges: drop merge of f%d (tested? %d) into f%d (tested? %d) mergetype %d, dist %2.2g, angle %4.4g. Merge independent sets of coplanar merges\n", facet1->id, facet1->visible, facet2->id, facet2->visible, mergetype, distance, angle)); continue; } } } trace3((qh ferr, 3047, "qh_all_merges: merge f%d and f%d type %d dist %2.2g angle %4.4g\n", facet1->id, facet2->id, mergetype, distance, angle)); if (mergetype == MRGtwisted) qh_merge_twisted(facet1, facet2); else qh_merge_nonconvex(facet1, facet2, mergetype); numnewmerges++; numdegenredun += qh_merge_degenredundant(); wasmerge= True; if (mergetype == MRGconcave) numconcave++; else if (mergetype == MRGconcavecoplanar) numconcavecoplanar++; else if (mergetype == MRGtwisted) numtwisted++; else if (mergetype == MRGcoplanar || mergetype == MRGanglecoplanar) numcoplanar++; else { qh_fprintf(qh ferr, 6394, "qhull internal error (qh_all_merges): expecting concave, coplanar, or twisted merge. Got merge f%d f%d v%d mergetype %d\n", getid_(facet1), getid_(facet2), getid_(vertex), mergetype); qh_errexit2(qh_ERRqhull, facet1, facet2); } } /* while qh_setdellast */ if (qh POSTmerging && qh hull_dim <= qh_DIMreduceBuild && numnewmerges > qh_MAXnewmerges) { numnewmerges= 0; wasmerge= othermerge= False; qh_reducevertices(); /* otherwise large post merges too slow */ } qh_getmergeset(qh newfacet_list); /* qh.facet_mergeset */ } /* while facet_mergeset or degen_mergeset */ if (qh VERTEXneighbors) { /* at least one merge */ isreduce= False; if (qh POSTmerging && qh hull_dim >= 4) { isreduce= True; }else if (qh POSTmerging || !qh MERGEexact) { if ((wasmerge || othermerge) && qh hull_dim > 2 && qh hull_dim <= qh_DIMreduceBuild) isreduce= True; } if (isreduce) { wasmerge= othermerge= False; if (qh_reducevertices()) { qh_getmergeset(qh newfacet_list); /* facet_mergeset */ continue; } } } if (vneighbors && qh_test_vneighbors(/* qh.newfacet_list */)) continue; break; } /* while (True) */ if (wasmerge || othermerge) { trace3((qh ferr, 3033, "qh_all_merges: skip qh_reducevertices due to post-merging, no qh.VERTEXneighbors (%d), or hull_dim %d ==2 or >%d\n", qh VERTEXneighbors, qh hull_dim, qh_DIMreduceBuild)) FORALLnew_facets { newfacet->newmerge= False; } } if (qh CHECKfrequently && !qh MERGEexact) { qh old_randomdist= qh RANDOMdist; qh RANDOMdist= False; qh_checkconvex(qh newfacet_list, qh_ALGORITHMfault); /* qh_checkconnect(); [this is slow and it changes the facet order] */ qh RANDOMdist= qh old_randomdist; } trace1((qh ferr, 1009, "qh_all_merges: merged %d coplanar %d concave %d concavecoplanar %d twisted facets and %d degen or redundant facets.\n", numcoplanar, numconcave, numconcavecoplanar, numtwisted, numdegenredun)); if (qh IStracing >= 4 && qh num_facets < 500) qh_printlists(); } /* all_merges */ /*--------------------------------- qh_all_vertexmerges( apexpointid, facet, &retryfacet ) merge vertices in qh.vertex_mergeset and subsequent merges returns: returns retryfacet for facet (if defined) updates qh.facet_list, qh.num_facets, qh.vertex_list, qh.num_vertices mergesets are empty if merges, resets facet lists notes: called from qh_qhull, qh_addpoint, and qh_buildcone_mergepinched vertex merges occur after facet merges and qh_resetlists design: while merges in vertex_mergeset (MRGvertices) merge a pair of pinched vertices update vertex neighbors merge non-convex and degenerate facets and check for ridges with duplicate vertices partition outside points of deleted, "visible" facets */ void qh_all_vertexmerges(int apexpointid, facetT *facet, facetT **retryfacet) { int numpoints; /* ignore count of partitioned points. Used by qh_addpoint for Zpbalance */ if (retryfacet) *retryfacet= facet; while (qh_setsize(qh vertex_mergeset) > 0) { trace1((qh ferr, 1057, "qh_all_vertexmerges: starting to merge %d vertex merges for apex p%d facet f%d\n", qh_setsize(qh vertex_mergeset), apexpointid, getid_(facet))); if (qh IStracing >= 4 && qh num_facets < 1000) qh_printlists(); qh_merge_pinchedvertices(apexpointid /* qh.vertex_mergeset, visible_list, newvertex_list, newfacet_list */); qh_update_vertexneighbors(); /* update neighbors of qh.newvertex_list from qh_newvertices for deleted facets on qh.visible_list */ /* test ridges and merge non-convex facets */ qh_getmergeset(qh newfacet_list); qh_all_merges(True, False); /* calls qh_reducevertices */ if (qh CHECKfrequently) qh_checkpolygon(qh facet_list); qh_partitionvisible(!qh_ALL, &numpoints /* qh.visible_list qh.del_vertices*/); if (retryfacet) *retryfacet= qh_getreplacement(*retryfacet); qh_deletevisible(/* qh.visible_list qh.del_vertices*/); qh_resetlists(False, qh_RESETvisible /* qh.visible_list newvertex_list qh.newfacet_list */); if (qh IStracing >= 4 && qh num_facets < 1000) { qh_printlists(); qh_checkpolygon(qh facet_list); } } } /* all_vertexmerges */ /*--------------------------------- qh_appendmergeset( facet, vertex, neighbor, mergetype, dist, angle ) appends an entry to qh.facet_mergeset or qh.degen_mergeset if 'dist' is unknown, set it to 0.0 if 'angle' is unknown, set it to 1.0 (coplanar) returns: merge appended to facet_mergeset or degen_mergeset sets ->degenerate or ->redundant if degen_mergeset notes: caller collects statistics and/or caller of qh_mergefacet see: qh_test_appendmerge() design: allocate merge entry if regular merge append to qh.facet_mergeset else if degenerate merge and qh.facet_mergeset is all degenerate append to qh.degen_mergeset else if degenerate merge prepend to qh.degen_mergeset (merged last) else if redundant merge append to qh.degen_mergeset */ void qh_appendmergeset(facetT *facet, facetT *neighbor, mergeType mergetype, coordT dist, realT angle) { mergeT *merge, *lastmerge; void **freelistp; /* used if !qh_NOmem by qh_memalloc_() */ const char *mergename; if ((facet->redundant && mergetype != MRGmirror) || neighbor->redundant) { trace3((qh ferr, 3051, "qh_appendmergeset: f%d is already redundant (%d) or f%d is already redundant (%d). Ignore merge f%d and f%d type %d\n", facet->id, facet->redundant, neighbor->id, neighbor->redundant, facet->id, neighbor->id, mergetype)); return; } if (facet->degenerate && mergetype == MRGdegen) { trace3((qh ferr, 3077, "qh_appendmergeset: f%d is already degenerate. Ignore merge f%d type %d (MRGdegen)\n", facet->id, facet->id, mergetype)); return; } if (!qh facet_mergeset || !qh degen_mergeset) { qh_fprintf(qh ferr, 6403, "qhull internal error (qh_appendmergeset): expecting temp set defined for qh.facet_mergeset (0x%x) and qh.degen_mergeset (0x%x). Got NULL\n", qh facet_mergeset, qh degen_mergeset); /* otherwise qh_setappend creates a new set that is not freed by qh_freebuild() */ qh_errexit(qh_ERRqhull, NULL, NULL); } if (neighbor->flipped && !facet->flipped) { if (mergetype != MRGdupridge) { qh_fprintf(qh ferr, 6355, "qhull internal error (qh_appendmergeset): except for MRGdupridge, cannot merge a non-flipped facet f%d into flipped f%d, mergetype %d, dist %4.4g\n", facet->id, neighbor->id, mergetype, dist); qh_errexit(qh_ERRqhull, NULL, NULL); }else { trace2((qh ferr, 2106, "qh_appendmergeset: dupridge will merge a non-flipped facet f%d into flipped f%d, dist %4.4g\n", facet->id, neighbor->id, dist)); } } qh_memalloc_((int)sizeof(mergeT), freelistp, merge, mergeT); merge->angle= angle; merge->distance= dist; merge->facet1= facet; merge->facet2= neighbor; merge->vertex1= NULL; merge->vertex2= NULL; merge->ridge1= NULL; merge->ridge2= NULL; merge->mergetype= mergetype; if(mergetype > 0 && mergetype < sizeof(mergetypes)/sizeof(char *)) mergename= mergetypes[mergetype]; else mergename= mergetypes[MRGnone]; if (mergetype < MRGdegen) qh_setappend(&(qh facet_mergeset), merge); else if (mergetype == MRGdegen) { facet->degenerate= True; if (!(lastmerge= (mergeT *)qh_setlast(qh degen_mergeset)) || lastmerge->mergetype == MRGdegen) qh_setappend(&(qh degen_mergeset), merge); else qh_setaddnth(&(qh degen_mergeset), 0, merge); /* merged last */ }else if (mergetype == MRGredundant) { facet->redundant= True; qh_setappend(&(qh degen_mergeset), merge); }else /* mergetype == MRGmirror */ { if (facet->redundant || neighbor->redundant) { qh_fprintf(qh ferr, 6092, "qhull internal error (qh_appendmergeset): facet f%d or f%d is already a mirrored facet (i.e., 'redundant')\n", facet->id, neighbor->id); qh_errexit2(qh_ERRqhull, facet, neighbor); } if (!qh_setequal(facet->vertices, neighbor->vertices)) { qh_fprintf(qh ferr, 6093, "qhull internal error (qh_appendmergeset): mirrored facets f%d and f%d do not have the same vertices\n", facet->id, neighbor->id); qh_errexit2(qh_ERRqhull, facet, neighbor); } facet->redundant= True; neighbor->redundant= True; qh_setappend(&(qh degen_mergeset), merge); } if (merge->mergetype >= MRGdegen) { trace3((qh ferr, 3044, "qh_appendmergeset: append merge f%d and f%d type %d (%s) to qh.degen_mergeset (size %d)\n", merge->facet1->id, merge->facet2->id, merge->mergetype, mergename, qh_setsize(qh degen_mergeset))); }else { trace3((qh ferr, 3027, "qh_appendmergeset: append merge f%d and f%d type %d (%s) dist %2.2g angle %4.4g to qh.facet_mergeset (size %d)\n", merge->facet1->id, merge->facet2->id, merge->mergetype, mergename, merge->distance, merge->angle, qh_setsize(qh facet_mergeset))); } } /* appendmergeset */ /*--------------------------------- qh_appendvertexmerge( vertex, vertex2, mergetype, distance, ridge1, ridge2 ) appends a vertex merge to qh.vertex_mergeset MRGsubridge includes two ridges (from MRGdupridge) MRGvertices includes two ridges notes: called by qh_getpinchedmerges for MRGsubridge called by qh_maybe_duplicateridge and qh_maybe_duplicateridges for MRGvertices only way to add a vertex merge to qh.vertex_mergeset checked by qh_next_vertexmerge */ void qh_appendvertexmerge(vertexT *vertex, vertexT *destination, mergeType mergetype, realT distance, ridgeT *ridge1, ridgeT *ridge2) { mergeT *merge; void **freelistp; /* used if !qh_NOmem by qh_memalloc_() */ const char *mergename; if (!qh vertex_mergeset) { qh_fprintf(qh ferr, 6387, "qhull internal error (qh_appendvertexmerge): expecting temp set defined for qh.vertex_mergeset (0x%x). Got NULL\n", qh vertex_mergeset); /* otherwise qh_setappend creates a new set that is not freed by qh_freebuild() */ qh_errexit(qh_ERRqhull, NULL, NULL); } qh_memalloc_((int)sizeof(mergeT), freelistp, merge, mergeT); merge->angle= qh_ANGLEnone; merge->distance= distance; merge->facet1= NULL; merge->facet2= NULL; merge->vertex1= vertex; merge->vertex2= destination; merge->ridge1= ridge1; merge->ridge2= ridge2; merge->mergetype= mergetype; if(mergetype > 0 && mergetype < sizeof(mergetypes)/sizeof(char *)) mergename= mergetypes[mergetype]; else mergename= mergetypes[MRGnone]; if (mergetype == MRGvertices) { if (!ridge1 || !ridge2 || ridge1 == ridge2) { qh_fprintf(qh ferr, 6106, "qhull internal error (qh_appendvertexmerge): expecting two distinct ridges for MRGvertices. Got r%d r%d\n", getid_(ridge1), getid_(ridge2)); qh_errexit(qh_ERRqhull, NULL, ridge1); } } qh_setappend(&(qh vertex_mergeset), merge); trace3((qh ferr, 3034, "qh_appendvertexmerge: append merge v%d into v%d r%d r%d dist %2.2g type %d (%s)\n", vertex->id, destination->id, getid_(ridge1), getid_(ridge2), distance, merge->mergetype, mergename)); } /* appendvertexmerge */ /*--------------------------------- qh_basevertices( samecycle ) return temporary set of base vertices for samecycle samecycle is first facet in the cycle assumes apex is SETfirst_( samecycle->vertices ) returns: vertices(settemp) all ->seen are cleared notes: uses qh_vertex_visit; design: for each facet in samecycle for each unseen vertex in facet->vertices append to result */ setT *qh_basevertices(facetT *samecycle) { facetT *same; vertexT *apex, *vertex, **vertexp; setT *vertices= qh_settemp(qh TEMPsize); apex= SETfirstt_(samecycle->vertices, vertexT); apex->visitid= ++qh vertex_visit; FORALLsame_cycle_(samecycle) { if (same->mergeridge) continue; FOREACHvertex_(same->vertices) { if (vertex->visitid != qh vertex_visit) { qh_setappend(&vertices, vertex); vertex->visitid= qh vertex_visit; vertex->seen= False; } } } trace4((qh ferr, 4019, "qh_basevertices: found %d vertices\n", qh_setsize(vertices))); return vertices; } /* basevertices */ /*--------------------------------- qh_check_dupridge( facet1, dist1, facet2, dist2 ) Check dupridge between facet1 and facet2 for wide merge dist1 is the maximum distance of facet1's vertices to facet2 dist2 is the maximum distance of facet2's vertices to facet1 returns Level 1 log of the dupridge with the minimum distance between vertices Throws error if the merge will increase the maximum facet width by qh_WIDEduplicate (100x) notes: only called from qh_forcedmerges */ void qh_check_dupridge(facetT *facet1, realT dist1, facetT *facet2, realT dist2) { vertexT *vertex, **vertexp, *vertexA, **vertexAp; realT dist, innerplane, mergedist, outerplane, prevdist, ratio, vertexratio; realT minvertex= REALmax; mergedist= fmin_(dist1, dist2); qh_outerinner(NULL, &outerplane, &innerplane); /* ratio from qh_printsummary */ FOREACHvertex_(facet1->vertices) { /* The dupridge is between facet1 and facet2, so either facet can be tested */ FOREACHvertexA_(facet1->vertices) { if (vertex > vertexA){ /* Test each pair once */ dist= qh_pointdist(vertex->point, vertexA->point, qh hull_dim); minimize_(minvertex, dist); /* Not quite correct. A facet may have a dupridge and another pair of nearly adjacent vertices. */ } } } prevdist= fmax_(outerplane, innerplane); maximize_(prevdist, qh ONEmerge + qh DISTround); maximize_(prevdist, qh MINoutside + qh DISTround); ratio= mergedist/prevdist; vertexratio= minvertex/prevdist; trace0((qh ferr, 16, "qh_check_dupridge: dupridge between f%d and f%d (vertex dist %2.2g), dist %2.2g, reverse dist %2.2g, ratio %2.2g while processing p%d\n", facet1->id, facet2->id, minvertex, dist1, dist2, ratio, qh furthest_id)); if (ratio > qh_WIDEduplicate) { qh_fprintf(qh ferr, 6271, "qhull topology error (qh_check_dupridge): wide merge (%.1fx wider) due to dupridge between f%d and f%d (vertex dist %2.2g), merge dist %2.2g, while processing p%d\n- Allow error with option 'Q12'\n", ratio, facet1->id, facet2->id, minvertex, mergedist, qh furthest_id); if (vertexratio < qh_WIDEpinched) qh_fprintf(qh ferr, 8145, "- Experimental option merge-pinched-vertices ('Q14') may avoid this error. It merges nearly adjacent vertices.\n"); if (qh DELAUNAY) qh_fprintf(qh ferr, 8145, "- A bounding box for the input sites may alleviate this error.\n"); if (!qh ALLOWwide) qh_errexit2(qh_ERRwide, facet1, facet2); } } /* check_dupridge */ /*--------------------------------- qh_checkconnect( ) check that new facets are connected new facets are on qh.newfacet_list notes: this is slow and it changes the order of the facets uses qh.visit_id design: move first new facet to end of qh.facet_list for all newly appended facets append unvisited neighbors to end of qh.facet_list for all new facets report error if unvisited */ void qh_checkconnect(void /* qh.newfacet_list */) { facetT *facet, *newfacet, *errfacet= NULL, *neighbor, **neighborp; facet= qh newfacet_list; qh_removefacet(facet); qh_appendfacet(facet); facet->visitid= ++qh visit_id; FORALLfacet_(facet) { FOREACHneighbor_(facet) { if (neighbor->visitid != qh visit_id) { qh_removefacet(neighbor); qh_appendfacet(neighbor); neighbor->visitid= qh visit_id; } } } FORALLnew_facets { if (newfacet->visitid == qh visit_id) break; qh_fprintf(qh ferr, 6094, "qhull internal error (qh_checkconnect): f%d is not attached to the new facets\n", newfacet->id); errfacet= newfacet; } if (errfacet) qh_errexit(qh_ERRqhull, errfacet, NULL); } /* checkconnect */ /*--------------------------------- qh_checkdelfacet( facet, mergeset ) check that mergeset does not reference facet */ void qh_checkdelfacet(facetT *facet, setT *mergeset) { mergeT *merge, **mergep; FOREACHmerge_(mergeset) { if (merge->facet1 == facet || merge->facet2 == facet) { qh_fprintf(qh ferr, 6390, "qhull internal error (qh_checkdelfacet): cannot delete f%d. It is referenced by merge f%d f%d mergetype %d\n", facet->id, merge->facet1->id, getid_(merge->facet2), merge->mergetype); qh_errexit2(qh_ERRqhull, merge->facet1, merge->facet2); } } } /* checkdelfacet */ /*--------------------------------- qh_checkdelridge( ) check that qh_delridge_merge is not needed for deleted ridges notes: called from qh_mergecycle, qh_makenewfacets, qh_attachnewfacets errors if qh.vertex_mergeset is non-empty errors if any visible or new facet has a ridge with r.nonconvex set assumes that vertex.delfacet is not needed */ void qh_checkdelridge(void /* qh.visible_facets, vertex_mergeset */) { facetT *newfacet, *visible; ridgeT *ridge, **ridgep; if (!SETempty_(qh vertex_mergeset)) { qh_fprintf(qh ferr, 6382, "qhull internal error (qh_checkdelridge): expecting empty qh.vertex_mergeset in order to avoid calling qh_delridge_merge. Got %d merges\n", qh_setsize(qh vertex_mergeset)); qh_errexit(qh_ERRqhull, NULL, NULL); } FORALLnew_facets { FOREACHridge_(newfacet->ridges) { if (ridge->nonconvex) { qh_fprintf(qh ferr, 6313, "qhull internal error (qh_checkdelridge): unexpected 'nonconvex' flag for ridge r%d in newfacet f%d. Otherwise need to call qh_delridge_merge\n", ridge->id, newfacet->id); qh_errexit(qh_ERRqhull, newfacet, ridge); } } } FORALLvisible_facets { FOREACHridge_(visible->ridges) { if (ridge->nonconvex) { qh_fprintf(qh ferr, 6385, "qhull internal error (qh_checkdelridge): unexpected 'nonconvex' flag for ridge r%d in visible facet f%d. Otherwise need to call qh_delridge_merge\n", ridge->id, visible->id); qh_errexit(qh_ERRqhull, visible, ridge); } } } } /* checkdelridge */ /*--------------------------------- qh_checkzero( testall ) check that facets are clearly convex for qh.DISTround with qh.MERGEexact if testall, test all facets for qh.MERGEexact post-merging else test qh.newfacet_list if qh.MERGEexact, allows coplanar ridges skips convexity test while qh.ZEROall_ok returns: True if all facets !flipped, !dupridge, normal if all horizon facets are simplicial if all vertices are clearly below neighbor if all opposite vertices of horizon are below clears qh.ZEROall_ok if any problems or coplanar facets notes: called by qh_premerge (qh.CHECKzero, 'C-0') and qh_qhull ('Qx') uses qh.vertex_visit horizon facets may define multiple new facets design: for all facets in qh.newfacet_list or qh.facet_list check for flagged faults (flipped, etc.) for all facets in qh.newfacet_list or qh.facet_list for each neighbor of facet skip horizon facets for qh.newfacet_list test the opposite vertex if qh.newfacet_list test the other vertices in the facet's horizon facet */ boolT qh_checkzero(boolT testall) { facetT *facet, *neighbor; facetT *horizon, *facetlist; int neighbor_i, neighbor_n; vertexT *vertex, **vertexp; realT dist; if (testall) facetlist= qh facet_list; else { facetlist= qh newfacet_list; FORALLfacet_(facetlist) { horizon= SETfirstt_(facet->neighbors, facetT); if (!horizon->simplicial) goto LABELproblem; if (facet->flipped || facet->dupridge || !facet->normal) goto LABELproblem; } if (qh MERGEexact && qh ZEROall_ok) { trace2((qh ferr, 2011, "qh_checkzero: skip convexity check until first pre-merge\n")); return True; } } FORALLfacet_(facetlist) { qh vertex_visit++; horizon= NULL; FOREACHneighbor_i_(facet) { if (!neighbor_i && !testall) { horizon= neighbor; continue; /* horizon facet tested in qh_findhorizon */ } vertex= SETelemt_(facet->vertices, neighbor_i, vertexT); vertex->visitid= qh vertex_visit; zzinc_(Zdistzero); qh_distplane(vertex->point, neighbor, &dist); if (dist >= -2 * qh DISTround) { /* need 2x for qh_distround and 'Rn' for qh_checkconvex, same as qh.premerge_centrum */ qh ZEROall_ok= False; if (!qh MERGEexact || testall || dist > qh DISTround) goto LABELnonconvex; } } if (!testall && horizon) { FOREACHvertex_(horizon->vertices) { if (vertex->visitid != qh vertex_visit) { zzinc_(Zdistzero); qh_distplane(vertex->point, facet, &dist); if (dist >= -2 * qh DISTround) { qh ZEROall_ok= False; if (!qh MERGEexact || dist > qh DISTround) goto LABELnonconvexhorizon; } break; } } } } trace2((qh ferr, 2012, "qh_checkzero: testall %d, facets are %s\n", testall, (qh MERGEexact && !testall) ? "not concave, flipped, or dupridge" : "clearly convex")); return True; LABELproblem: qh ZEROall_ok= False; trace2((qh ferr, 2013, "qh_checkzero: qh_premerge is needed. New facet f%d or its horizon f%d is non-simplicial, flipped, dupridge, or mergehorizon\n", facet->id, horizon->id)); return False; LABELnonconvex: trace2((qh ferr, 2014, "qh_checkzero: facet f%d and f%d are not clearly convex. v%d dist %.2g\n", facet->id, neighbor->id, vertex->id, dist)); return False; LABELnonconvexhorizon: trace2((qh ferr, 2060, "qh_checkzero: facet f%d and horizon f%d are not clearly convex. v%d dist %.2g\n", facet->id, horizon->id, vertex->id, dist)); return False; } /* checkzero */ /*--------------------------------- qh_compare_anglemerge( mergeA, mergeB ) used by qsort() to order qh.facet_mergeset by mergetype and angle (qh.ANGLEmerge, 'Q1') lower numbered mergetypes done first (MRGcoplanar before MRGconcave) notes: qh_all_merges processes qh.facet_mergeset by qh_setdellast [mar'19] evaluated various options with eg/q_benchmark and merging of pinched vertices (Q14) */ int qh_compare_anglemerge(const void *p1, const void *p2) { const mergeT *a= *((mergeT *const*)p1), *b= *((mergeT *const*)p2); if (a->mergetype != b->mergetype) return (a->mergetype < b->mergetype ? 1 : -1); /* select MRGcoplanar (1) before MRGconcave (3) */ else return (a->angle > b->angle ? 1 : -1); /* select coplanar merge (1.0) before sharp merge (-0.5) */ } /* compare_anglemerge */ /*--------------------------------- qh_compare_facetmerge( mergeA, mergeB ) used by qsort() to order merges by mergetype, first merge, first lower numbered mergetypes done first (MRGcoplanar before MRGconcave) if same merge type, flat merges are first notes: qh_all_merges processes qh.facet_mergeset by qh_setdellast [mar'19] evaluated various options with eg/q_benchmark and merging of pinched vertices (Q14) */ int qh_compare_facetmerge(const void *p1, const void *p2) { const mergeT *a= *((mergeT *const*)p1), *b= *((mergeT *const*)p2); if (a->mergetype != b->mergetype) return (a->mergetype < b->mergetype ? 1 : -1); /* select MRGcoplanar (1) before MRGconcave (3) */ else if (a->mergetype == MRGanglecoplanar) return (a->angle > b->angle ? 1 : -1); /* if MRGanglecoplanar, select coplanar merge (1.0) before sharp merge (-0.5) */ else return (a->distance < b->distance ? 1 : -1); /* select flat (0.0) merge before wide (1e-10) merge */ } /* compare_facetmerge */ /*--------------------------------- qh_comparevisit( vertexA, vertexB ) used by qsort() to order vertices by their visitid notes: only called by qh_find_newvertex */ int qh_comparevisit(const void *p1, const void *p2) { const vertexT *a= *((vertexT *const*)p1), *b= *((vertexT *const*)p2); if (a->visitid > b->visitid) return 1; return -1; } /* comparevisit */ /*--------------------------------- qh_copynonconvex( atridge ) set non-convex flag on other ridges (if any) between same neighbors notes: may be faster if use smaller ridge set design: for each ridge of atridge's top facet if ridge shares the same neighbor set nonconvex flag */ void qh_copynonconvex(ridgeT *atridge) { facetT *facet, *otherfacet; ridgeT *ridge, **ridgep; facet= atridge->top; otherfacet= atridge->bottom; atridge->nonconvex= False; FOREACHridge_(facet->ridges) { if (otherfacet == ridge->top || otherfacet == ridge->bottom) { if (ridge != atridge) { ridge->nonconvex= True; trace4((qh ferr, 4020, "qh_copynonconvex: moved nonconvex flag from r%d to r%d between f%d and f%d\n", atridge->id, ridge->id, facet->id, otherfacet->id)); break; } } } } /* copynonconvex */ /*--------------------------------- qh_degen_redundant_facet( facet ) check for a degenerate (too few neighbors) or redundant (subset of vertices) facet notes: called at end of qh_mergefacet, qh_renamevertex, and qh_reducevertices bumps vertex_visit called if a facet was redundant but no longer is (qh_merge_degenredundant) qh_appendmergeset() only appends first reference to facet (i.e., redundant) see: qh_test_redundant_neighbors, qh_maydropneighbor design: test for redundant neighbor test for degenerate facet */ void qh_degen_redundant_facet(facetT *facet) { vertexT *vertex, **vertexp; facetT *neighbor, **neighborp; trace3((qh ferr, 3028, "qh_degen_redundant_facet: test facet f%d for degen/redundant\n", facet->id)); if (facet->flipped) { trace2((qh ferr, 3074, "qh_degen_redundant_facet: f%d is flipped, will merge later\n", facet->id)); return; } FOREACHneighbor_(facet) { if (neighbor->flipped) /* disallow merge of non-flipped into flipped, neighbor will be merged later */ continue; if (neighbor->visible) { qh_fprintf(qh ferr, 6357, "qhull internal error (qh_degen_redundant_facet): facet f%d has deleted neighbor f%d (qh.visible_list)\n", facet->id, neighbor->id); qh_errexit2(qh_ERRqhull, facet, neighbor); } qh vertex_visit++; FOREACHvertex_(neighbor->vertices) vertex->visitid= qh vertex_visit; FOREACHvertex_(facet->vertices) { if (vertex->visitid != qh vertex_visit) break; } if (!vertex) { trace2((qh ferr, 2015, "qh_degen_redundant_facet: f%d is contained in f%d. merge\n", facet->id, neighbor->id)); qh_appendmergeset(facet, neighbor, MRGredundant, 0.0, 1.0); return; } } if (qh_setsize(facet->neighbors) < qh hull_dim) { qh_appendmergeset(facet, facet, MRGdegen, 0.0, 1.0); trace2((qh ferr, 2016, "qh_degen_redundant_facet: f%d is degenerate.\n", facet->id)); } } /* degen_redundant_facet */ /*--------------------------------- qh_delridge_merge( ridge ) delete ridge due to a merge notes: only called by merge.c (qh_mergeridges, qh_renameridgevertex) ridges also freed in qh_freeqhull and qh_mergecycle_ridges design: if needed, moves ridge.nonconvex to another ridge sets vertex.delridge for qh_reducevertices deletes ridge from qh.vertex_mergeset deletes ridge from its neighboring facets frees up its memory */ void qh_delridge_merge(ridgeT *ridge) { vertexT *vertex, **vertexp; mergeT *merge; int merge_i, merge_n; trace3((qh ferr, 3036, "qh_delridge_merge: delete ridge r%d between f%d and f%d\n", ridge->id, ridge->top->id, ridge->bottom->id)); if (ridge->nonconvex) qh_copynonconvex(ridge); FOREACHvertex_(ridge->vertices) vertex->delridge= True; FOREACHmerge_i_(qh vertex_mergeset) { if (merge->ridge1 == ridge || merge->ridge2 == ridge) { trace3((qh ferr, 3029, "qh_delridge_merge: drop merge of v%d into v%d (dist %2.2g r%d r%d) due to deleted, duplicated ridge r%d\n", merge->vertex1->id, merge->vertex2->id, merge->distance, merge->ridge1->id, merge->ridge2->id, ridge->id)); if (merge->ridge1 == ridge) merge->ridge2->mergevertex= False; else merge->ridge1->mergevertex= False; qh_setdelnth(qh vertex_mergeset, merge_i); merge_i--; merge_n--; /* next merge after deleted */ } } qh_setdel(ridge->top->ridges, ridge); qh_setdel(ridge->bottom->ridges, ridge); qh_delridge(ridge); } /* delridge_merge */ /*--------------------------------- qh_drop_mergevertex( merge ) clear mergevertex flags for ridges of a vertex merge */ void qh_drop_mergevertex(mergeT *merge) { if (merge->mergetype == MRGvertices) { merge->ridge1->mergevertex= False; merge->ridge1->mergevertex2= True; merge->ridge2->mergevertex= False; merge->ridge2->mergevertex2= True; trace3((qh ferr, 3032, "qh_drop_mergevertex: unset mergevertex for r%d and r%d due to dropped vertex merge v%d to v%d. Sets mergevertex2\n", merge->ridge1->id, merge->ridge2->id, merge->vertex1->id, merge->vertex2->id)); } } /* drop_mergevertex */ /*--------------------------------- qh_find_newvertex( oldvertex, vertices, ridges ) locate new vertex for renaming old vertex vertices is a set of possible new vertices vertices sorted by number of deleted ridges returns: newvertex or NULL each ridge includes both newvertex and oldvertex vertices without oldvertex sorted by number of deleted ridges qh.vertex_visit updated sets v.seen notes: called by qh_redundant_vertex due to vertex->delridge and qh_rename_sharedvertex sets vertex->visitid to 0..setsize() for vertices new vertex is in one of the ridges renaming will not cause a duplicate ridge renaming will minimize the number of deleted ridges newvertex may not be adjacent in the dual (though unlikely) design: for each vertex in vertices set vertex->visitid to number of ridges remove unvisited vertices set qh.vertex_visit above all possible values sort vertices by number of ridges (minimize ridges that need renaming add each ridge to qh.hash_table for each vertex in vertices find the first vertex that would not cause a duplicate ridge after a rename */ vertexT *qh_find_newvertex(vertexT *oldvertex, setT *vertices, setT *ridges) { vertexT *vertex, **vertexp; setT *newridges; ridgeT *ridge, **ridgep; int size, hashsize; int hash; unsigned int maxvisit; #ifndef qh_NOtrace if (qh IStracing >= 4) { qh_fprintf(qh ferr, 8063, "qh_find_newvertex: find new vertex for v%d from ", oldvertex->id); FOREACHvertex_(vertices) qh_fprintf(qh ferr, 8064, "v%d ", vertex->id); FOREACHridge_(ridges) qh_fprintf(qh ferr, 8065, "r%d ", ridge->id); qh_fprintf(qh ferr, 8066, "\n"); } #endif FOREACHridge_(ridges) { FOREACHvertex_(ridge->vertices) vertex->seen= False; } FOREACHvertex_(vertices) { vertex->visitid= 0; /* v.visitid will be number of ridges */ vertex->seen= True; } FOREACHridge_(ridges) { FOREACHvertex_(ridge->vertices) { if (vertex->seen) vertex->visitid++; } } FOREACHvertex_(vertices) { if (!vertex->visitid) { qh_setdelnth(vertices, SETindex_(vertices,vertex)); vertexp--; /* repeat since deleted this vertex */ } } maxvisit= (unsigned int)qh_setsize(ridges); maximize_(qh vertex_visit, maxvisit); if (!qh_setsize(vertices)) { trace4((qh ferr, 4023, "qh_find_newvertex: vertices not in ridges for v%d\n", oldvertex->id)); return NULL; } qsort(SETaddr_(vertices, vertexT), (size_t)qh_setsize(vertices), sizeof(vertexT *), qh_comparevisit); /* can now use qh vertex_visit */ if (qh PRINTstatistics) { size= qh_setsize(vertices); zinc_(Zintersect); zadd_(Zintersecttot, size); zmax_(Zintersectmax, size); } hashsize= qh_newhashtable(qh_setsize(ridges)); FOREACHridge_(ridges) qh_hashridge(qh hash_table, hashsize, ridge, oldvertex); FOREACHvertex_(vertices) { newridges= qh_vertexridges(vertex, !qh_ALL); FOREACHridge_(newridges) { if (qh_hashridge_find(qh hash_table, hashsize, ridge, vertex, oldvertex, &hash)) { zinc_(Zvertexridge); break; } } qh_settempfree(&newridges); if (!ridge) break; /* found a rename */ } if (vertex) { /* counted in qh_renamevertex */ trace2((qh ferr, 2020, "qh_find_newvertex: found v%d for old v%d from %d vertices and %d ridges.\n", vertex->id, oldvertex->id, qh_setsize(vertices), qh_setsize(ridges))); }else { zinc_(Zfindfail); trace0((qh ferr, 14, "qh_find_newvertex: no vertex for renaming v%d (all duplicated ridges) during p%d\n", oldvertex->id, qh furthest_id)); } qh_setfree(&qh hash_table); return vertex; } /* find_newvertex */ /*--------------------------------- qh_findbest_pinchedvertex( merge, apex, nearestp, distp ) Determine the best pinched vertex to rename as its nearest neighboring vertex Renaming will remove a duplicate MRGdupridge in newfacet_list returns: pinched vertex (either apex or subridge), nearest vertex (subridge or neighbor vertex), and the distance between them notes: only called by qh_getpinchedmerges assumes qh.VERTEXneighbors see qh_findbest_ridgevertex design: if the facets have the same vertices return the nearest vertex pair else the subridge is the intersection of the two new facets minus the apex the subridge consists of qh.hull_dim-2 horizon vertices the subridge is also a matched ridge for the new facets (its duplicate) determine the nearest vertex to the apex determine the nearest pair of subridge vertices for each vertex in the subridge determine the nearest neighbor vertex (not in the subridge) */ vertexT *qh_findbest_pinchedvertex(mergeT *merge, vertexT *apex, vertexT **nearestp, coordT *distp /* qh.newfacet_list */) { vertexT *vertex, **vertexp, *vertexA, **vertexAp; vertexT *bestvertex= NULL, *bestpinched= NULL; setT *subridge, *maybepinched; coordT dist, bestdist= REALmax; coordT pincheddist= (qh ONEmerge+qh DISTround)*qh_RATIOpinchedsubridge; if (!merge->facet1->simplicial || !merge->facet2->simplicial) { qh_fprintf(qh ferr, 6351, "qhull internal error (qh_findbest_pinchedvertex): expecting merge of adjacent, simplicial new facets. f%d or f%d is not simplicial\n", merge->facet1->id, merge->facet2->id); qh_errexit2(qh_ERRqhull, merge->facet1, merge->facet2); } subridge= qh_vertexintersect_new(merge->facet1->vertices, merge->facet2->vertices); /* new setT. No error_exit() */ if (qh_setsize(subridge) == qh hull_dim) { /* duplicate vertices */ bestdist= qh_vertex_bestdist2(subridge, &bestvertex, &bestpinched); if(bestvertex == apex) { bestvertex= bestpinched; bestpinched= apex; } }else { qh_setdel(subridge, apex); if (qh_setsize(subridge) != qh hull_dim - 2) { qh_fprintf(qh ferr, 6409, "qhull internal error (qh_findbest_pinchedvertex): expecting subridge of qh.hull_dim-2 vertices for the intersection of new facets f%d and f%d minus their apex. Got %d vertices\n", merge->facet1->id, merge->facet2->id, qh_setsize(subridge)); qh_errexit2(qh_ERRqhull, merge->facet1, merge->facet2); } FOREACHvertex_(subridge) { dist= qh_pointdist(vertex->point, apex->point, qh hull_dim); if (dist < bestdist) { bestpinched= apex; bestvertex= vertex; bestdist= dist; } } if (bestdist > pincheddist) { FOREACHvertex_(subridge) { FOREACHvertexA_(subridge) { if (vertexA->id > vertex->id) { /* once per vertex pair, do not compare addresses */ dist= qh_pointdist(vertexA->point, vertex->point, qh hull_dim); if (dist < bestdist) { bestpinched= vertexA; bestvertex= vertex; bestdist= dist; } } } } } if (bestdist > pincheddist) { FOREACHvertexA_(subridge) { maybepinched= qh_neighbor_vertices(vertexA, subridge); /* subridge and apex tested above */ FOREACHvertex_(maybepinched) { dist= qh_pointdist(vertex->point, vertexA->point, qh hull_dim); if (dist < bestdist) { bestvertex= vertex; bestpinched= vertexA; bestdist= dist; } } qh_settempfree(&maybepinched); } } } *distp= bestdist; qh_setfree(&subridge); /* qh_err_exit not called since allocated */ if (!bestvertex) { /* should never happen if qh.hull_dim > 2 */ qh_fprintf(qh ferr, 6274, "qhull internal error (qh_findbest_pinchedvertex): did not find best vertex for subridge of dupridge between f%d and f%d, while processing p%d\n", merge->facet1->id, merge->facet2->id, qh furthest_id); qh_errexit2(qh_ERRqhull, merge->facet1, merge->facet2); } *nearestp= bestvertex; trace2((qh ferr, 2061, "qh_findbest_pinchedvertex: best pinched p%d(v%d) and vertex p%d(v%d) are closest (%2.2g) for duplicate subridge between f%d and f%d\n", qh_pointid(bestpinched->point), bestpinched->id, qh_pointid(bestvertex->point), bestvertex->id, bestdist, merge->facet1->id, merge->facet2->id)); return bestpinched; } /* findbest_pinchedvertex */ /*--------------------------------- qh_findbest_ridgevertex( ridge, pinchedp, distp ) Determine the best vertex/pinched-vertex to merge for ridges with the same vertices returns: vertex, pinched vertex, and the distance between them notes: assumes qh.hull_dim>=3 see qh_findbest_pinchedvertex */ vertexT *qh_findbest_ridgevertex(ridgeT *ridge, vertexT **pinchedp, coordT *distp) { vertexT *bestvertex; *distp= qh_vertex_bestdist2(ridge->vertices, &bestvertex, pinchedp); trace4((qh ferr, 4069, "qh_findbest_ridgevertex: best pinched p%d(v%d) and vertex p%d(v%d) are closest (%2.2g) for duplicated ridge r%d (same vertices) between f%d and f%d\n", qh_pointid((*pinchedp)->point), (*pinchedp)->id, qh_pointid(bestvertex->point), bestvertex->id, *distp, ridge->id, ridge->top->id, ridge->bottom->id)); return bestvertex; } /* findbest_ridgevertex */ /*--------------------------------- qh_findbest_test( testcentrum, facet, neighbor, &bestfacet, &dist, &mindist, &maxdist ) test neighbor of facet for qh_findbestneighbor() if testcentrum, tests centrum (assumes it is defined) else tests vertices initially *bestfacet==NULL and *dist==REALmax returns: if a better facet (i.e., vertices/centrum of facet closer to neighbor) updates bestfacet, dist, mindist, and maxdist notes: called by qh_findbestneighbor ignores pairs of flipped facets, unless that's all there is */ void qh_findbest_test(boolT testcentrum, facetT *facet, facetT *neighbor, facetT **bestfacet, realT *distp, realT *mindistp, realT *maxdistp) { realT dist, mindist, maxdist; if (facet->flipped && neighbor->flipped && *bestfacet && !(*bestfacet)->flipped) return; /* do not merge flipped into flipped facets */ if (testcentrum) { zzinc_(Zbestdist); qh_distplane(facet->center, neighbor, &dist); dist *= qh hull_dim; /* estimate furthest vertex */ if (dist < 0) { maxdist= 0; mindist= dist; dist= -dist; }else { mindist= 0; maxdist= dist; } }else dist= qh_getdistance(facet, neighbor, &mindist, &maxdist); if (dist < *distp) { *bestfacet= neighbor; *mindistp= mindist; *maxdistp= maxdist; *distp= dist; } } /* findbest_test */ /*--------------------------------- qh_findbestneighbor( facet, dist, mindist, maxdist ) finds best neighbor (least dist) of a facet for merging returns: returns min and max distances and their max absolute value notes: error if qh_ASvoronoi avoids merging old into new assumes ridge->nonconvex only set on one ridge between a pair of facets could use an early out predicate but not worth it design: if a large facet will test centrum else will test vertices if a large facet test nonconvex neighbors for best merge else test all neighbors for the best merge if testing centrum get distance information */ facetT *qh_findbestneighbor(facetT *facet, realT *distp, realT *mindistp, realT *maxdistp) { facetT *neighbor, **neighborp, *bestfacet= NULL; ridgeT *ridge, **ridgep; boolT nonconvex= True, testcentrum= False; int size= qh_setsize(facet->vertices); if(qh CENTERtype==qh_ASvoronoi){ qh_fprintf(qh ferr, 6272, "qhull internal error: cannot call qh_findbestneighor for f%d while qh.CENTERtype is qh_ASvoronoi\n", facet->id); qh_errexit(qh_ERRqhull, facet, NULL); } *distp= REALmax; if (size > qh_BESTcentrum2 * qh hull_dim + qh_BESTcentrum) { testcentrum= True; zinc_(Zbestcentrum); if (!facet->center) facet->center= qh_getcentrum(facet); } if (size > qh hull_dim + qh_BESTnonconvex) { FOREACHridge_(facet->ridges) { if (ridge->nonconvex) { neighbor= otherfacet_(ridge, facet); qh_findbest_test(testcentrum, facet, neighbor, &bestfacet, distp, mindistp, maxdistp); } } } if (!bestfacet) { nonconvex= False; FOREACHneighbor_(facet) qh_findbest_test(testcentrum, facet, neighbor, &bestfacet, distp, mindistp, maxdistp); } if (!bestfacet) { qh_fprintf(qh ferr, 6095, "qhull internal error (qh_findbestneighbor): no neighbors for f%d\n", facet->id); qh_errexit(qh_ERRqhull, facet, NULL); } if (testcentrum) qh_getdistance(facet, bestfacet, mindistp, maxdistp); trace3((qh ferr, 3002, "qh_findbestneighbor: f%d is best neighbor for f%d testcentrum? %d nonconvex? %d dist %2.2g min %2.2g max %2.2g\n", bestfacet->id, facet->id, testcentrum, nonconvex, *distp, *mindistp, *maxdistp)); return(bestfacet); } /* findbestneighbor */ /*--------------------------------- qh_flippedmerges( facetlist, wasmerge ) merge flipped facets into best neighbor assumes qh.facet_mergeset at top of temporary stack returns: no flipped facets on facetlist sets wasmerge if merge occurred degen/redundant merges passed through notes: othermerges not needed since qh.facet_mergeset is empty before & after keep it in case of change design: append flipped facets to qh.facetmergeset for each flipped merge find best neighbor merge facet into neighbor merge degenerate and redundant facets remove flipped merges from qh.facet_mergeset */ void qh_flippedmerges(facetT *facetlist, boolT *wasmerge) { facetT *facet, *neighbor, *facet1; realT dist, mindist, maxdist; mergeT *merge, **mergep; setT *othermerges; int nummerge= 0, numdegen= 0; trace4((qh ferr, 4024, "qh_flippedmerges: begin\n")); FORALLfacet_(facetlist) { if (facet->flipped && !facet->visible) qh_appendmergeset(facet, facet, MRGflip, 0.0, 1.0); } othermerges= qh_settemppop(); if(othermerges != qh facet_mergeset) { qh_fprintf(qh ferr, 6392, "qhull internal error (qh_flippedmerges): facet_mergeset (%d merges) not at top of tempstack (%d merges)\n", qh_setsize(qh facet_mergeset), qh_setsize(othermerges)); qh_errexit(qh_ERRqhull, NULL, NULL); } qh facet_mergeset= qh_settemp(qh TEMPsize); qh_settemppush(othermerges); FOREACHmerge_(othermerges) { facet1= merge->facet1; if (merge->mergetype != MRGflip || facet1->visible) continue; if (qh TRACEmerge-1 == zzval_(Ztotmerge)) qhmem.IStracing= qh IStracing= qh TRACElevel; neighbor= qh_findbestneighbor(facet1, &dist, &mindist, &maxdist); trace0((qh ferr, 15, "qh_flippedmerges: merge flipped f%d into f%d dist %2.2g during p%d\n", facet1->id, neighbor->id, dist, qh furthest_id)); qh_mergefacet(facet1, neighbor, merge->mergetype, &mindist, &maxdist, !qh_MERGEapex); nummerge++; if (qh PRINTstatistics) { zinc_(Zflipped); wadd_(Wflippedtot, dist); wmax_(Wflippedmax, dist); } } FOREACHmerge_(othermerges) { if (merge->facet1->visible || merge->facet2->visible) qh_memfree(merge, (int)sizeof(mergeT)); /* invalidates merge and othermerges */ else qh_setappend(&qh facet_mergeset, merge); } qh_settempfree(&othermerges); numdegen += qh_merge_degenredundant(); /* somewhat better here than after each flipped merge -- qtest.sh 10 '500 C1,2e-13 D4' 'd Qbb' */ if (nummerge) *wasmerge= True; trace1((qh ferr, 1010, "qh_flippedmerges: merged %d flipped and %d degenredundant facets into a good neighbor\n", nummerge, numdegen)); } /* flippedmerges */ /*--------------------------------- qh_forcedmerges( wasmerge ) merge dupridges calls qh_check_dupridge to report an error on wide merges assumes qh_settemppop is qh.facet_mergeset returns: removes all dupridges on facet_mergeset wasmerge set if merge qh.facet_mergeset may include non-forced merges(none for now) qh.degen_mergeset includes degen/redun merges notes: called by qh_premerge dupridges occur when the horizon is pinched, i.e. a subridge occurs in more than two horizon ridges. could rename vertices that pinch the horizon assumes qh_merge_degenredundant() has not be called othermerges isn't needed since facet_mergeset is empty afterwards keep it in case of change design: for each dupridge find current facets by chasing f.replace links check for wide merge due to dupridge determine best direction for facet merge one facet into the other remove dupridges from qh.facet_mergeset */ void qh_forcedmerges(boolT *wasmerge) { facetT *facet1, *facet2, *merging, *merged, *newfacet; mergeT *merge, **mergep; realT dist, mindist, maxdist, dist2, mindist2, maxdist2; setT *othermerges; int nummerge=0, numflip=0, numdegen= 0; boolT wasdupridge= False; if (qh TRACEmerge-1 == zzval_(Ztotmerge)) qhmem.IStracing= qh IStracing= qh TRACElevel; trace3((qh ferr, 3054, "qh_forcedmerges: merge dupridges\n")); othermerges= qh_settemppop(); /* was facet_mergeset */ if (qh facet_mergeset != othermerges ) { qh_fprintf(qh ferr, 6279, "qhull internal error (qh_forcedmerges): qh_settemppop (size %d) is not qh facet_mergeset (size %d)\n", qh_setsize(othermerges), qh_setsize(qh facet_mergeset)); qh_errexit(qh_ERRqhull, NULL, NULL); } qh facet_mergeset= qh_settemp(qh TEMPsize); qh_settemppush(othermerges); FOREACHmerge_(othermerges) { if (merge->mergetype != MRGdupridge) continue; wasdupridge= True; if (qh TRACEmerge-1 == zzval_(Ztotmerge)) qhmem.IStracing= qh IStracing= qh TRACElevel; facet1= qh_getreplacement(merge->facet1); /* must exist, no qh_merge_degenredunant */ facet2= qh_getreplacement(merge->facet2); /* previously merged facet, if any */ if (facet1 == facet2) continue; if (!qh_setin(facet2->neighbors, facet1)) { qh_fprintf(qh ferr, 6096, "qhull internal error (qh_forcedmerges): f%d and f%d had a dupridge but as f%d and f%d they are no longer neighbors\n", merge->facet1->id, merge->facet2->id, facet1->id, facet2->id); qh_errexit2(qh_ERRqhull, facet1, facet2); } dist= qh_getdistance(facet1, facet2, &mindist, &maxdist); dist2= qh_getdistance(facet2, facet1, &mindist2, &maxdist2); qh_check_dupridge(facet1, dist, facet2, dist2); if (dist < dist2) { if (facet2->flipped && !facet1->flipped && dist2 < qh_WIDEdupridge*(qh ONEmerge+qh DISTround)) { /* prefer merge of flipped facet */ merging= facet2; merged= facet1; dist= dist2; mindist= mindist2; maxdist= maxdist2; }else { merging= facet1; merged= facet2; } }else { if (facet1->flipped && !facet2->flipped && dist < qh_WIDEdupridge*(qh ONEmerge+qh DISTround)) { /* prefer merge of flipped facet */ merging= facet1; merged= facet2; }else { merging= facet2; merged= facet1; dist= dist2; mindist= mindist2; maxdist= maxdist2; } } qh_mergefacet(merging, merged, merge->mergetype, &mindist, &maxdist, !qh_MERGEapex); numdegen += qh_merge_degenredundant(); /* better here than at end -- qtest.sh 10 '500 C1,2e-13 D4' 'd Qbb' */ if (facet1->flipped) { zinc_(Zmergeflipdup); numflip++; }else nummerge++; if (qh PRINTstatistics) { zinc_(Zduplicate); wadd_(Wduplicatetot, dist); wmax_(Wduplicatemax, dist); } } FOREACHmerge_(othermerges) { if (merge->mergetype == MRGdupridge) qh_memfree(merge, (int)sizeof(mergeT)); /* invalidates merge and othermerges */ else qh_setappend(&qh facet_mergeset, merge); } qh_settempfree(&othermerges); if (wasdupridge) { FORALLnew_facets { if (newfacet->dupridge) { newfacet->dupridge= False; newfacet->mergeridge= False; newfacet->mergeridge2= False; if (qh_setsize(newfacet->neighbors) < qh hull_dim) { /* not tested for MRGdupridge */ qh_appendmergeset(newfacet, newfacet, MRGdegen, 0.0, 1.0); trace2((qh ferr, 2107, "qh_forcedmerges: dupridge f%d is degenerate with fewer than %d neighbors\n", newfacet->id, qh hull_dim)); } } } numdegen += qh_merge_degenredundant(); } if (nummerge || numflip) { *wasmerge= True; trace1((qh ferr, 1011, "qh_forcedmerges: merged %d facets, %d flipped facets, and %d degenredundant facets across dupridges\n", nummerge, numflip, numdegen)); } } /* forcedmerges */ /*--------------------------------- qh_freemergesets( ) free the merge sets notes: matches qh_initmergesets */ void qh_freemergesets(void) { if (!qh facet_mergeset || !qh degen_mergeset || !qh vertex_mergeset) { qh_fprintf(qh ferr, 6388, "qhull internal error (qh_freemergesets): expecting mergesets. Got a NULL mergeset, qh.facet_mergeset (0x%x), qh.degen_mergeset (0x%x), qh.vertex_mergeset (0x%x)\n", qh facet_mergeset, qh degen_mergeset, qh vertex_mergeset); qh_errexit(qh_ERRqhull, NULL, NULL); } if (!SETempty_(qh facet_mergeset) || !SETempty_(qh degen_mergeset) || !SETempty_(qh vertex_mergeset)) { qh_fprintf(qh ferr, 6389, "qhull internal error (qh_freemergesets): expecting empty mergesets. Got qh.facet_mergeset (%d merges), qh.degen_mergeset (%d merges), qh.vertex_mergeset (%d merges)\n", qh_setsize(qh facet_mergeset), qh_setsize(qh degen_mergeset), qh_setsize(qh vertex_mergeset)); qh_errexit(qh_ERRqhull, NULL, NULL); } qh_settempfree(&qh facet_mergeset); qh_settempfree(&qh vertex_mergeset); qh_settempfree(&qh degen_mergeset); } /* freemergesets */ /*--------------------------------- qh_getmergeset( facetlist ) determines nonconvex facets on facetlist tests !tested ridges and nonconvex ridges of !tested facets returns: returns sorted qh.facet_mergeset of facet-neighbor pairs to be merged all ridges tested notes: facetlist is qh.facet_newlist, use qh_getmergeset_initial for all facets assumes no nonconvex ridges with both facets tested uses facet->tested/ridge->tested to prevent duplicate tests can not limit tests to modified ridges since the centrum changed uses qh.visit_id design: for each facet on facetlist for each ridge of facet if untested ridge test ridge for convexity if non-convex append ridge to qh.facet_mergeset sort qh.facet_mergeset by mergetype and angle or distance */ void qh_getmergeset(facetT *facetlist) { facetT *facet, *neighbor, **neighborp; ridgeT *ridge, **ridgep; int nummerges; boolT simplicial; nummerges= qh_setsize(qh facet_mergeset); trace4((qh ferr, 4026, "qh_getmergeset: started.\n")); qh visit_id++; FORALLfacet_(facetlist) { if (facet->tested) continue; facet->visitid= qh visit_id; FOREACHneighbor_(facet) neighbor->seen= False; /* facet must be non-simplicial due to merge to qh.facet_newlist */ FOREACHridge_(facet->ridges) { if (ridge->tested && !ridge->nonconvex) continue; /* if r.tested & r.nonconvex, need to retest and append merge */ neighbor= otherfacet_(ridge, facet); if (neighbor->seen) { /* another ridge for this facet-neighbor pair was already tested in this loop */ ridge->tested= True; ridge->nonconvex= False; /* only one ridge is marked nonconvex per facet-neighbor pair */ }else if (neighbor->visitid != qh visit_id) { neighbor->seen= True; ridge->nonconvex= False; simplicial= False; if (ridge->simplicialbot && ridge->simplicialtop) simplicial= True; if (qh_test_appendmerge(facet, neighbor, simplicial)) ridge->nonconvex= True; ridge->tested= True; } } facet->tested= True; } nummerges= qh_setsize(qh facet_mergeset); if (qh ANGLEmerge) qsort(SETaddr_(qh facet_mergeset, mergeT), (size_t)nummerges, sizeof(mergeT *), qh_compare_anglemerge); else qsort(SETaddr_(qh facet_mergeset, mergeT), (size_t)nummerges, sizeof(mergeT *), qh_compare_facetmerge); nummerges += qh_setsize(qh degen_mergeset); if (qh POSTmerging) { zadd_(Zmergesettot2, nummerges); }else { zadd_(Zmergesettot, nummerges); zmax_(Zmergesetmax, nummerges); } trace2((qh ferr, 2021, "qh_getmergeset: %d merges found\n", nummerges)); } /* getmergeset */ /*--------------------------------- qh_getmergeset_initial( facetlist ) determine initial qh.facet_mergeset for facets tests all facet/neighbor pairs on facetlist returns: sorted qh.facet_mergeset with nonconvex ridges sets facet->tested, ridge->tested, and ridge->nonconvex notes: uses visit_id, assumes ridge->nonconvex is False see qh_getmergeset design: for each facet on facetlist for each untested neighbor of facet test facet and neighbor for convexity if non-convex append merge to qh.facet_mergeset mark one of the ridges as nonconvex sort qh.facet_mergeset by mergetype and angle or distance */ void qh_getmergeset_initial(facetT *facetlist) { facetT *facet, *neighbor, **neighborp; ridgeT *ridge, **ridgep; int nummerges; boolT simplicial; qh visit_id++; FORALLfacet_(facetlist) { facet->visitid= qh visit_id; FOREACHneighbor_(facet) { if (neighbor->visitid != qh visit_id) { simplicial= False; /* ignores r.simplicialtop/simplicialbot. Need to test horizon facets */ if (facet->simplicial && neighbor->simplicial) simplicial= True; if (qh_test_appendmerge(facet, neighbor, simplicial)) { FOREACHridge_(neighbor->ridges) { if (facet == otherfacet_(ridge, neighbor)) { ridge->nonconvex= True; break; /* only one ridge is marked nonconvex */ } } } } } facet->tested= True; FOREACHridge_(facet->ridges) ridge->tested= True; } nummerges= qh_setsize(qh facet_mergeset); if (qh ANGLEmerge) qsort(SETaddr_(qh facet_mergeset, mergeT), (size_t)nummerges, sizeof(mergeT *), qh_compare_anglemerge); else qsort(SETaddr_(qh facet_mergeset, mergeT), (size_t)nummerges, sizeof(mergeT *), qh_compare_facetmerge); nummerges += qh_setsize(qh degen_mergeset); if (qh POSTmerging) { zadd_(Zmergeinittot2, nummerges); }else { zadd_(Zmergeinittot, nummerges); zmax_(Zmergeinitmax, nummerges); } trace2((qh ferr, 2022, "qh_getmergeset_initial: %d merges found\n", nummerges)); } /* getmergeset_initial */ /*--------------------------------- qh_getpinchedmerges( apex, maxdist, iscoplanar ) get pinched merges for dupridges in qh.facet_mergeset qh.NEWtentative==True qh.newfacet_list with apex qh.horizon_list is attached to qh.visible_list instead of qh.newfacet_list maxdist for vertex-facet of a dupridge qh.facet_mergeset is empty qh.vertex_mergeset is a temporary set returns: False if nearest vertex would increase facet width by more than maxdist or qh_WIDEpinched True and iscoplanar, if the pinched vertex is the apex (i.e., make the apex a coplanar point) True and !iscoplanar, if should merge a pinched vertex of a dupridge qh.vertex_mergeset contains one or more MRGsubridge with a pinched vertex and a nearby, neighboring vertex qh.facet_mergeset is empty notes: called by qh_buildcone_mergepinched hull_dim >= 3 a pinched vertex is in a dupridge and the horizon selects the pinched vertex that is closest to its neighbor design: for each dupridge determine the best pinched vertex to be merged into a neighboring vertex if merging the pinched vertex would produce a wide merge (qh_WIDEpinched) ignore pinched vertex with a warning, and use qh_merge_degenredundant instead else append the pinched vertex to vertex_mergeset for merging */ boolT qh_getpinchedmerges(vertexT *apex, coordT maxdupdist, boolT *iscoplanar /* qh.newfacet_list, qh.vertex_mergeset */) { mergeT *merge, **mergep, *bestmerge= NULL; vertexT *nearest, *pinched, *bestvertex= NULL, *bestpinched= NULL; boolT result; coordT dist, prevdist, bestdist= REALmax/(qh_RATIOcoplanarapex+1.0); /* allow *3.0 */ realT ratio; trace2((qh ferr, 2062, "qh_getpinchedmerges: try to merge pinched vertices for dupridges in new facets with apex p%d(v%d) max dupdist %2.2g\n", qh_pointid(apex->point), apex->id, maxdupdist)); *iscoplanar= False; prevdist= fmax_(qh ONEmerge + qh DISTround, qh MINoutside + qh DISTround); maximize_(prevdist, qh max_outside); maximize_(prevdist, -qh min_vertex); qh_mark_dupridges(qh newfacet_list, !qh_ALL); /* qh.facet_mergeset, creates ridges */ /* qh_mark_dupridges is called a second time in qh_premerge */ FOREACHmerge_(qh facet_mergeset) { /* read-only */ if (merge->mergetype != MRGdupridge) { qh_fprintf(qh ferr, 6393, "qhull internal error (qh_getpinchedmerges): expecting MRGdupridge from qh_mark_dupridges. Got merge f%d f%d type %d\n", getid_(merge->facet1), getid_(merge->facet2), merge->mergetype); qh_errexit(qh_ERRqhull, NULL, NULL); } /* dist is distance between vertices */ pinched= qh_findbest_pinchedvertex(merge, apex, &nearest, &dist /* qh.newfacet_list */); if (pinched == apex && dist < qh_RATIOcoplanarapex*bestdist) { /* prefer coplanar apex since it always works */ bestdist= dist/qh_RATIOcoplanarapex; bestmerge= merge; bestpinched= pinched; bestvertex= nearest; }else if (dist < bestdist) { bestdist= dist; bestmerge= merge; bestpinched= pinched; bestvertex= nearest; } } result= False; if (bestmerge && bestdist < maxdupdist) { ratio= bestdist / prevdist; if (ratio > qh_WIDEpinched) { if (bestmerge->facet1->mergehorizon || bestmerge->facet2->mergehorizon) { /* e.g., rbox 175 C3,2e-13 t1539182828 | qhull d */ trace1((qh ferr, 1051, "qh_getpinchedmerges: dupridge (MRGdupridge) of coplanar horizon would produce a wide merge (%.0fx) due to pinched vertices v%d and v%d (dist %2.2g) for f%d and f%d. qh_mergecycle_all will merge one or both facets\n", ratio, bestpinched->id, bestvertex->id, bestdist, bestmerge->facet1->id, bestmerge->facet2->id)); }else { qh_fprintf(qh ferr, 7081, "qhull precision warning (qh_getpinchedmerges): pinched vertices v%d and v%d (dist %2.2g, %.0fx) would produce a wide merge for f%d and f%d. Will merge dupridge instead\n", bestpinched->id, bestvertex->id, bestdist, ratio, bestmerge->facet1->id, bestmerge->facet2->id); } }else { if (bestpinched == apex) { trace2((qh ferr, 2063, "qh_getpinchedmerges: will make the apex a coplanar point. apex p%d(v%d) is the nearest vertex to v%d on dupridge. Dist %2.2g\n", qh_pointid(apex->point), apex->id, bestvertex->id, bestdist*qh_RATIOcoplanarapex)); qh coplanar_apex= apex->point; *iscoplanar= True; result= True; }else if (qh_setin(bestmerge->facet1->vertices, bestpinched) != qh_setin(bestmerge->facet2->vertices, bestpinched)) { /* pinched in one facet but not the other facet */ trace2((qh ferr, 2064, "qh_getpinchedmerges: will merge new facets to resolve dupridge between f%d and f%d with pinched v%d and v%d\n", bestmerge->facet1->id, bestmerge->facet2->id, bestpinched->id, bestvertex->id)); qh_appendvertexmerge(bestpinched, bestvertex, MRGsubridge, bestdist, NULL, NULL); result= True; }else { trace2((qh ferr, 2065, "qh_getpinchedmerges: will merge pinched v%d into v%d to resolve dupridge between f%d and f%d\n", bestpinched->id, bestvertex->id, bestmerge->facet1->id, bestmerge->facet2->id)); qh_appendvertexmerge(bestpinched, bestvertex, MRGsubridge, bestdist, NULL, NULL); result= True; } } } /* delete MRGdupridge, qh_mark_dupridges is called a second time in qh_premerge */ while ((merge= (mergeT *)qh_setdellast(qh facet_mergeset))) qh_memfree(merge, (int)sizeof(mergeT)); return result; }/* getpinchedmerges */ /*--------------------------------- qh_hasmerge( mergeset, mergetype, facetA, facetB ) True if mergeset has mergetype for facetA and facetB */ boolT qh_hasmerge(setT *mergeset, mergeType type, facetT *facetA, facetT *facetB) { mergeT *merge, **mergep; FOREACHmerge_(mergeset) { if (merge->mergetype == type) { if (merge->facet1 == facetA && merge->facet2 == facetB) return True; if (merge->facet1 == facetB && merge->facet2 == facetA) return True; } } return False; }/* hasmerge */ /*--------------------------------- qh_hashridge( hashtable, hashsize, ridge, oldvertex ) add ridge to hashtable without oldvertex notes: assumes hashtable is large enough design: determine hash value for ridge without oldvertex find next empty slot for ridge */ void qh_hashridge(setT *hashtable, int hashsize, ridgeT *ridge, vertexT *oldvertex) { int hash; ridgeT *ridgeA; hash= qh_gethash(hashsize, ridge->vertices, qh hull_dim-1, 0, oldvertex); while (True) { if (!(ridgeA= SETelemt_(hashtable, hash, ridgeT))) { SETelem_(hashtable, hash)= ridge; break; }else if (ridgeA == ridge) break; if (++hash == hashsize) hash= 0; } } /* hashridge */ /*--------------------------------- qh_hashridge_find( hashtable, hashsize, ridge, vertex, oldvertex, hashslot ) returns matching ridge without oldvertex in hashtable for ridge without vertex if oldvertex is NULL matches with any one skip returns: matching ridge or NULL if no match, if ridge already in table hashslot= -1 else hashslot= next NULL index notes: assumes hashtable is large enough can't match ridge to itself design: get hash value for ridge without vertex for each hashslot return match if ridge matches ridgeA without oldvertex */ ridgeT *qh_hashridge_find(setT *hashtable, int hashsize, ridgeT *ridge, vertexT *vertex, vertexT *oldvertex, int *hashslot) { int hash; ridgeT *ridgeA; *hashslot= 0; zinc_(Zhashridge); hash= qh_gethash(hashsize, ridge->vertices, qh hull_dim-1, 0, vertex); while ((ridgeA= SETelemt_(hashtable, hash, ridgeT))) { if (ridgeA == ridge) *hashslot= -1; else { zinc_(Zhashridgetest); if (qh_setequal_except(ridge->vertices, vertex, ridgeA->vertices, oldvertex)) return ridgeA; } if (++hash == hashsize) hash= 0; } if (!*hashslot) *hashslot= hash; return NULL; } /* hashridge_find */ /*--------------------------------- qh_initmergesets( ) initialize the merge sets if 'all', include qh.degen_mergeset notes: matches qh_freemergesets */ void qh_initmergesets(void /* qh.facet_mergeset,degen_mergeset,vertex_mergeset */) { if (qh facet_mergeset || qh degen_mergeset || qh vertex_mergeset) { qh_fprintf(qh ferr, 6386, "qhull internal error (qh_initmergesets): expecting NULL mergesets. Got qh.facet_mergeset (0x%x), qh.degen_mergeset (0x%x), qh.vertex_mergeset (0x%x)\n", qh facet_mergeset, qh degen_mergeset, qh vertex_mergeset); qh_errexit(qh_ERRqhull, NULL, NULL); } qh degen_mergeset= qh_settemp(qh TEMPsize); qh vertex_mergeset= qh_settemp(qh TEMPsize); qh facet_mergeset= qh_settemp(qh TEMPsize); /* last temporary set for qh_forcedmerges */ } /* initmergesets */ /*--------------------------------- qh_makeridges( facet ) creates explicit ridges between simplicial facets returns: facet with ridges and without qh_MERGEridge ->simplicial is False if facet was tested, new ridges are tested notes: allows qh_MERGEridge flag uses existing ridges duplicate neighbors ok if ridges already exist (qh_mergecycle_ridges) see: qh_mergecycle_ridges() qh_rename_adjacentvertex for qh_merge_pinchedvertices design: look for qh_MERGEridge neighbors mark neighbors that already have ridges for each unprocessed neighbor of facet create a ridge for neighbor and facet if any qh_MERGEridge neighbors delete qh_MERGEridge flags (previously processed by qh_mark_dupridges) */ void qh_makeridges(facetT *facet) { facetT *neighbor, **neighborp; ridgeT *ridge, **ridgep; int neighbor_i, neighbor_n; boolT toporient, mergeridge= False; if (!facet->simplicial) return; trace4((qh ferr, 4027, "qh_makeridges: make ridges for f%d\n", facet->id)); facet->simplicial= False; FOREACHneighbor_(facet) { if (neighbor == qh_MERGEridge) mergeridge= True; else neighbor->seen= False; } FOREACHridge_(facet->ridges) otherfacet_(ridge, facet)->seen= True; FOREACHneighbor_i_(facet) { if (neighbor == qh_MERGEridge) continue; /* fixed by qh_mark_dupridges */ else if (!neighbor->seen) { /* no current ridges */ ridge= qh_newridge(); ridge->vertices= qh_setnew_delnthsorted(facet->vertices, qh hull_dim, neighbor_i, 0); toporient= (boolT)(facet->toporient ^ (neighbor_i & 0x1)); if (toporient) { ridge->top= facet; ridge->bottom= neighbor; ridge->simplicialtop= True; ridge->simplicialbot= neighbor->simplicial; }else { ridge->top= neighbor; ridge->bottom= facet; ridge->simplicialtop= neighbor->simplicial; ridge->simplicialbot= True; } if (facet->tested && !mergeridge) ridge->tested= True; #if 0 /* this also works */ flip= (facet->toporient ^ neighbor->toporient)^(skip1 & 0x1) ^ (skip2 & 0x1); if (facet->toporient ^ (skip1 & 0x1) ^ flip) { ridge->top= neighbor; ridge->bottom= facet; ridge->simplicialtop= True; ridge->simplicialbot= neighbor->simplicial; }else { ridge->top= facet; ridge->bottom= neighbor; ridge->simplicialtop= neighbor->simplicial; ridge->simplicialbot= True; } #endif qh_setappend(&(facet->ridges), ridge); trace5((qh ferr, 5005, "makeridges: appended r%d to ridges for f%d. Next is ridges for neighbor f%d\n", ridge->id, facet->id, neighbor->id)); qh_setappend(&(neighbor->ridges), ridge); if (qh ridge_id == qh traceridge_id) qh traceridge= ridge; } } if (mergeridge) { while (qh_setdel(facet->neighbors, qh_MERGEridge)) ; /* delete each one */ } } /* makeridges */ /*--------------------------------- qh_mark_dupridges( facetlist, allmerges ) add duplicated ridges to qh.facet_mergeset facet-dupridge is true if it contains a subridge shared by more than one new facet for each such facet, one has a neighbor marked qh_MERGEridge allmerges is true if merging dupridges allmerges is false if merging pinched vertices followed by retry addpoint qh_mark_dupridges will be called again if pinched vertices not found returns: dupridges on qh.facet_mergeset (MRGdupridge) f.mergeridge and f.mergeridge2 set for facet f.mergeridge set for neighbor if allmerges is true make ridges for facets with dupridges as marked by qh_MERGEridge and both sides facet->dupridge removes qh_MERGEridge from neighbor sets notes: called by qh_premerge and qh_getpinchedmerges dupridges are due to duplicate subridges i.e. a subridge occurs in more than two horizon ridges. i.e., a ridge has more than two neighboring facets dupridges occur in at least two cases 1) a pinched horizon with nearly adjacent vertices -> merge the vertices (qh_getpinchedmerges) 2) more than one newfacet for a horizon face -> merge coplanar facets (qh_premerge) qh_matchdupridge previously identified the furthest apart pair of facets to retain they must have a matching subridge and the same orientation only way to set facet->mergeridge and mergeridge2 uses qh.visit_id design: for all facets on facetlist if facet contains a dupridge for each neighbor of facet if neighbor marked qh_MERGEridge (one side of the merge) set facet->mergeridge else if neighbor contains a dupridge and the back link is qh_MERGEridge append dupridge to qh.facet_mergeset exit if !allmerges for repeating qh_mark_dupridges later for each dupridge make ridge sets in preparation for merging remove qh_MERGEridge from neighbor set for each dupridge restore the missing neighbor from the neighbor set that was qh_MERGEridge add the missing ridge for this neighbor */ void qh_mark_dupridges(facetT *facetlist, boolT allmerges) { facetT *facet, *neighbor, **neighborp; int nummerge=0; mergeT *merge, **mergep; trace4((qh ferr, 4028, "qh_mark_dupridges: identify dupridges in facetlist f%d, allmerges? %d\n", facetlist->id, allmerges)); FORALLfacet_(facetlist) { /* not necessary for first call */ facet->mergeridge2= False; facet->mergeridge= False; } FORALLfacet_(facetlist) { if (facet->dupridge) { FOREACHneighbor_(facet) { if (neighbor == qh_MERGEridge) { facet->mergeridge= True; continue; } if (neighbor->dupridge) { if (!qh_setin(neighbor->neighbors, facet)) { /* i.e., it is qh_MERGEridge, neighbors are distinct */ qh_appendmergeset(facet, neighbor, MRGdupridge, 0.0, 1.0); facet->mergeridge2= True; facet->mergeridge= True; nummerge++; }else if (qh_setequal(facet->vertices, neighbor->vertices)) { /* neighbors are the same except for horizon and qh_MERGEridge, see QH7085 */ trace3((qh ferr, 3043, "qh_mark_dupridges): dupridge due to duplicate vertices for subridges f%d and f%d\n", facet->id, neighbor->id)); qh_appendmergeset(facet, neighbor, MRGdupridge, 0.0, 1.0); facet->mergeridge2= True; facet->mergeridge= True; nummerge++; break; /* same for all neighbors */ } } } } } if (!nummerge) return; if (!allmerges) { trace1((qh ferr, 1012, "qh_mark_dupridges: found %d duplicated ridges (MRGdupridge) for qh_getpinchedmerges\n", nummerge)); return; } trace1((qh ferr, 1048, "qh_mark_dupridges: found %d duplicated ridges (MRGdupridge) for qh_premerge. Prepare facets for merging\n", nummerge)); /* make ridges in preparation for merging */ FORALLfacet_(facetlist) { if (facet->mergeridge && !facet->mergeridge2) qh_makeridges(facet); } trace3((qh ferr, 3075, "qh_mark_dupridges: restore missing neighbors and ridges due to qh_MERGEridge\n")); FOREACHmerge_(qh facet_mergeset) { /* restore the missing neighbors */ if (merge->mergetype == MRGdupridge) { /* only between simplicial facets */ if (merge->facet2->mergeridge2 && qh_setin(merge->facet2->neighbors, merge->facet1)) { /* Due to duplicate or multiple subridges, e.g., ../eg/qtest.sh t712682 '200 s W1e-13 C1,1e-13 D5' 'd' merge->facet1: - neighboring facets: f27779 f59186 f59186 f59186 MERGEridge f59186 merge->facet2: - neighboring facets: f27779 f59100 f59100 f59100 f59100 f59100 or, ../eg/qtest.sh 100 '500 s W1e-13 C1,1e-13 D4' 'd' both facets will be degenerate after merge, consider for special case handling */ qh_fprintf(qh ferr, 6361, "qhull topological error (qh_mark_dupridges): multiple dupridges for f%d and f%d, including reverse\n", merge->facet1->id, merge->facet2->id); qh_errexit2(qh_ERRtopology, merge->facet1, merge->facet2); }else qh_setappend(&merge->facet2->neighbors, merge->facet1); qh_makeridges(merge->facet1); /* and the missing ridges */ } } } /* mark_dupridges */ /*--------------------------------- qh_maybe_duplicateridge( ridge ) add MRGvertices if neighboring facet has another ridge with the same vertices returns: adds rename requests to qh.vertex_mergeset notes: called by qh_renamevertex nop if 2-D expensive test Duplicate ridges may lead to new facets with same vertex set (QH7084), will try merging vertices same as qh_maybe_duplicateridges design: for the two neighbors if non-simplicial for each ridge with the same first and last vertices (max id and min id) if the remaining vertices are the same get the closest pair of vertices add to vertex_mergeset for merging */ void qh_maybe_duplicateridge(ridgeT *ridgeA) { ridgeT *ridge, **ridgep; vertexT *vertex, *pinched; facetT *neighbor; coordT dist; int i, k, last= qh hull_dim-2; if (qh hull_dim < 3 ) return; for (neighbor= ridgeA->top, i=0; i<2; neighbor= ridgeA->bottom, i++) { if (!neighbor->simplicial && neighbor->nummerge > 0) { /* skip degenerate neighbors with both new and old vertices that will be merged */ FOREACHridge_(neighbor->ridges) { if (ridge != ridgeA && SETfirst_(ridge->vertices) == SETfirst_(ridgeA->vertices)) { if (SETelem_(ridge->vertices, last) == SETelem_(ridgeA->vertices, last)) { for (k=1; kvertices, k) != SETelem_(ridgeA->vertices, k)) break; } if (k == last) { vertex= qh_findbest_ridgevertex(ridge, &pinched, &dist); trace2((qh ferr, 2069, "qh_maybe_duplicateridge: will merge v%d into v%d (dist %2.2g) due to duplicate ridges r%d/r%d with the same vertices. mergevertex set\n", pinched->id, vertex->id, dist, ridgeA->id, ridge->id, ridgeA->top->id, ridgeA->bottom->id, ridge->top->id, ridge->bottom->id)); qh_appendvertexmerge(pinched, vertex, MRGvertices, dist, ridgeA, ridge); ridge->mergevertex= True; /* disables check for duplicate vertices in qh_checkfacet */ ridgeA->mergevertex= True; } } } } } } } /* maybe_duplicateridge */ /*--------------------------------- qh_maybe_duplicateridges( facet ) if Q15, add MRGvertices if facet has ridges with the same vertices returns: adds rename requests to qh.vertex_mergeset notes: called at end of qh_mergefacet and qh_mergecycle_all only enabled if qh.CHECKduplicates ('Q15') and 3-D or more expensive test, not worth it same as qh_maybe_duplicateridge design: for all ridge pairs in facet if the same first and last vertices (max id and min id) if the remaining vertices are the same get the closest pair of vertices add to vertex_mergeset for merging */ void qh_maybe_duplicateridges(facetT *facet) { facetT *otherfacet; ridgeT *ridge, *ridge2; vertexT *vertex, *pinched; coordT dist; int ridge_i, ridge_n, i, k, last_v= qh hull_dim-2; if (qh hull_dim < 3 || !qh CHECKduplicates) return; FOREACHridge_i_(facet->ridges) { otherfacet= otherfacet_(ridge, facet); if (otherfacet->degenerate || otherfacet->redundant || otherfacet->dupridge || otherfacet->flipped) /* will merge */ continue; for (i=ridge_i+1; i < ridge_n; i++) { ridge2= SETelemt_(facet->ridges, i, ridgeT); otherfacet= otherfacet_(ridge2, facet); if (otherfacet->degenerate || otherfacet->redundant || otherfacet->dupridge || otherfacet->flipped) /* will merge */ continue; /* optimize qh_setequal(ridge->vertices, ridge2->vertices) */ if (SETelem_(ridge->vertices, last_v) == SETelem_(ridge2->vertices, last_v)) { /* SETfirst is likely to be the same */ if (SETfirst_(ridge->vertices) == SETfirst_(ridge2->vertices)) { for (k=1; k vertices, k) != SETelem_(ridge2->vertices, k)) break; } if (k == last_v) { vertex= qh_findbest_ridgevertex(ridge, &pinched, &dist); if (ridge->top == ridge2->bottom && ridge->bottom == ridge2->top) { /* proof that ridges may have opposite orientation */ trace2((qh ferr, 2088, "qh_maybe_duplicateridges: will merge v%d into v%d (dist %2.2g) due to opposite oriented ridges r%d/r%d for f%d and f%d\n", pinched->id, vertex->id, dist, ridge->id, ridge2->id, ridge->top->id, ridge->bottom->id)); }else { trace2((qh ferr, 2083, "qh_maybe_duplicateridges: will merge v%d into v%d (dist %2.2g) due to duplicate ridges with the same vertices r%d/r%d in merged facet f%d\n", pinched->id, vertex->id, dist, ridge->id, ridge2->id, facet->id)); } qh_appendvertexmerge(pinched, vertex, MRGvertices, dist, ridge, ridge2); ridge->mergevertex= True; /* disables check for duplicate vertices in qh_checkfacet */ ridge2->mergevertex= True; } } } } } } /* maybe_duplicateridges */ /*--------------------------------- qh_maydropneighbor( facet ) drop neighbor relationship if ridge was deleted between a non-simplicial facet and its neighbors returns: for deleted ridges ridges made for simplicial neighbors neighbor sets updated appends degenerate facets to qh.facet_mergeset notes: called by qh_renamevertex assumes neighbors do not include qh_MERGEridge (qh_makeridges) won't cause redundant facets since vertex inclusion is the same may drop vertex and neighbor if no ridge uses qh.visit_id design: visit all neighbors with ridges for each unvisited neighbor of facet delete neighbor and facet from the non-simplicial neighbor sets if neighbor becomes degenerate append neighbor to qh.degen_mergeset if facet is degenerate append facet to qh.degen_mergeset */ void qh_maydropneighbor(facetT *facet) { ridgeT *ridge, **ridgep; facetT *neighbor, **neighborp; qh visit_id++; trace4((qh ferr, 4029, "qh_maydropneighbor: test f%d for no ridges to a neighbor\n", facet->id)); if (facet->simplicial) { qh_fprintf(qh ferr, 6278, "qhull internal error (qh_maydropneighbor): not valid for simplicial f%d while adding furthest p%d\n", facet->id, qh furthest_id); qh_errexit(qh_ERRqhull, facet, NULL); } FOREACHridge_(facet->ridges) { ridge->top->visitid= qh visit_id; ridge->bottom->visitid= qh visit_id; } FOREACHneighbor_(facet) { if (neighbor->visible) { qh_fprintf(qh ferr, 6358, "qhull internal error (qh_maydropneighbor): facet f%d has deleted neighbor f%d (qh.visible_list)\n", facet->id, neighbor->id); qh_errexit2(qh_ERRqhull, facet, neighbor); } if (neighbor->visitid != qh visit_id) { trace2((qh ferr, 2104, "qh_maydropneighbor: facets f%d and f%d are no longer neighbors while adding furthest p%d\n", facet->id, neighbor->id, qh furthest_id)); if (neighbor->simplicial) { qh_fprintf(qh ferr, 6280, "qhull internal error (qh_maydropneighbor): not valid for simplicial neighbor f%d of f%d while adding furthest p%d\n", neighbor->id, facet->id, qh furthest_id); qh_errexit2(qh_ERRqhull, neighbor, facet); } zinc_(Zdropneighbor); qh_setdel(neighbor->neighbors, facet); if (qh_setsize(neighbor->neighbors) < qh hull_dim) { zinc_(Zdropdegen); qh_appendmergeset(neighbor, neighbor, MRGdegen, 0.0, qh_ANGLEnone); trace2((qh ferr, 2023, "qh_maydropneighbors: f%d is degenerate.\n", neighbor->id)); } qh_setdel(facet->neighbors, neighbor); neighborp--; /* repeat, deleted a neighbor */ } } if (qh_setsize(facet->neighbors) < qh hull_dim) { zinc_(Zdropdegen); qh_appendmergeset(facet, facet, MRGdegen, 0.0, qh_ANGLEnone); trace2((qh ferr, 2024, "qh_maydropneighbors: f%d is degenerate.\n", facet->id)); } } /* maydropneighbor */ /*--------------------------------- qh_merge_degenredundant( ) merge all degenerate and redundant facets qh.degen_mergeset contains merges from qh_test_degen_neighbors, qh_test_redundant_neighbors, and qh_degen_redundant_facet returns: number of merges performed resets facet->degenerate/redundant if deleted (visible) facet has no neighbors sets ->f.replace to NULL notes: redundant merges happen before degenerate ones merging and renaming vertices can result in degen/redundant facets check for coplanar and convex neighbors afterwards design: for each merge on qh.degen_mergeset if redundant merge if non-redundant facet merged into redundant facet recheck facet for redundancy else merge redundant facet into other facet */ int qh_merge_degenredundant(void) { int size; mergeT *merge; facetT *bestneighbor, *facet1, *facet2, *facet3; realT dist, mindist, maxdist; vertexT *vertex, **vertexp; int nummerges= 0; mergeType mergetype; setT *mergedfacets; trace2((qh ferr, 2095, "qh_merge_degenredundant: merge %d degenerate, redundant, and mirror facets\n", qh_setsize(qh degen_mergeset))); mergedfacets= qh_settemp(qh TEMPsize); while ((merge= (mergeT *)qh_setdellast(qh degen_mergeset))) { facet1= merge->facet1; facet2= merge->facet2; mergetype= merge->mergetype; qh_memfree(merge, (int)sizeof(mergeT)); /* 'merge' is invalidated */ if (facet1->visible) continue; facet1->degenerate= False; facet1->redundant= False; if (qh TRACEmerge-1 == zzval_(Ztotmerge)) qhmem.IStracing= qh IStracing= qh TRACElevel; if (mergetype == MRGredundant) { zinc_(Zredundant); facet3= qh_getreplacement(facet2); /* the same facet if !facet2.visible */ if (!facet3) { qh_fprintf(qh ferr, 6097, "qhull internal error (qh_merge_degenredunant): f%d is redundant but visible f%d has no replacement\n", facet1->id, getid_(facet2)); qh_errexit2(qh_ERRqhull, facet1, facet2); } qh_setunique(&mergedfacets, facet3); if (facet1 == facet3) { continue; } trace2((qh ferr, 2025, "qh_merge_degenredundant: merge redundant f%d into f%d (arg f%d)\n", facet1->id, facet3->id, facet2->id)); qh_mergefacet(facet1, facet3, mergetype, NULL, NULL, !qh_MERGEapex); /* merge distance is already accounted for */ nummerges++; }else { /* mergetype == MRGdegen or MRGmirror, other merges may have fixed */ if (!(size= qh_setsize(facet1->neighbors))) { zinc_(Zdelfacetdup); trace2((qh ferr, 2026, "qh_merge_degenredundant: facet f%d has no neighbors. Deleted\n", facet1->id)); qh_willdelete(facet1, NULL); FOREACHvertex_(facet1->vertices) { qh_setdel(vertex->neighbors, facet1); if (!SETfirst_(vertex->neighbors)) { zinc_(Zdegenvertex); trace2((qh ferr, 2027, "qh_merge_degenredundant: deleted v%d because f%d has no neighbors\n", vertex->id, facet1->id)); vertex->deleted= True; qh_setappend(&qh del_vertices, vertex); } } nummerges++; }else if (size < qh hull_dim) { bestneighbor= qh_findbestneighbor(facet1, &dist, &mindist, &maxdist); trace2((qh ferr, 2028, "qh_merge_degenredundant: facet f%d has %d neighbors, merge into f%d dist %2.2g\n", facet1->id, size, bestneighbor->id, dist)); qh_mergefacet(facet1, bestneighbor, mergetype, &mindist, &maxdist, !qh_MERGEapex); nummerges++; if (qh PRINTstatistics) { zinc_(Zdegen); wadd_(Wdegentot, dist); wmax_(Wdegenmax, dist); } } /* else, another merge fixed the degeneracy and redundancy tested */ } } qh_settempfree(&mergedfacets); return nummerges; } /* merge_degenredundant */ /*--------------------------------- qh_merge_nonconvex( facet1, facet2, mergetype ) remove non-convex ridge between facet1 into facet2 mergetype gives why the facet's are non-convex returns: merges one of the facets into the best neighbor notes: mergetype is MRGcoplanar..MRGconvex design: if one of the facets is a new facet prefer merging new facet into old facet find best neighbors for both facets merge the nearest facet into its best neighbor update the statistics */ void qh_merge_nonconvex(facetT *facet1, facetT *facet2, mergeType mergetype) { facetT *bestfacet, *bestneighbor, *neighbor, *merging, *merged; realT dist, dist2, mindist, mindist2, maxdist, maxdist2; if (mergetype < MRGcoplanar || mergetype > MRGconcavecoplanar) { qh_fprintf(qh ferr, 6398, "qhull internal error (qh_merge_nonconvex): expecting mergetype MRGcoplanar..MRGconcavecoplanar. Got merge f%d and f%d type %d\n", facet1->id, facet2->id, mergetype); qh_errexit2(qh_ERRqhull, facet1, facet2); } if (qh TRACEmerge-1 == zzval_(Ztotmerge)) qhmem.IStracing= qh IStracing= qh TRACElevel; trace3((qh ferr, 3003, "qh_merge_nonconvex: merge #%d for f%d and f%d type %d\n", zzval_(Ztotmerge) + 1, facet1->id, facet2->id, mergetype)); /* concave or coplanar */ if (!facet1->newfacet) { bestfacet= facet2; /* avoid merging old facet if new is ok */ facet2= facet1; facet1= bestfacet; }else bestfacet= facet1; bestneighbor= qh_findbestneighbor(bestfacet, &dist, &mindist, &maxdist); neighbor= qh_findbestneighbor(facet2, &dist2, &mindist2, &maxdist2); if (dist < dist2) { merging= bestfacet; merged= bestneighbor; }else if (qh AVOIDold && !facet2->newfacet && ((mindist >= -qh MAXcoplanar && maxdist <= qh max_outside) || dist * 1.5 < dist2)) { zinc_(Zavoidold); wadd_(Wavoidoldtot, dist); wmax_(Wavoidoldmax, dist); trace2((qh ferr, 2029, "qh_merge_nonconvex: avoid merging old facet f%d dist %2.2g. Use f%d dist %2.2g instead\n", facet2->id, dist2, facet1->id, dist2)); merging= bestfacet; merged= bestneighbor; }else { merging= facet2; merged= neighbor; dist= dist2; mindist= mindist2; maxdist= maxdist2; } qh_mergefacet(merging, merged, mergetype, &mindist, &maxdist, !qh_MERGEapex); /* caller merges qh_degenredundant */ if (qh PRINTstatistics) { if (mergetype == MRGanglecoplanar) { zinc_(Zacoplanar); wadd_(Wacoplanartot, dist); wmax_(Wacoplanarmax, dist); }else if (mergetype == MRGconcave) { zinc_(Zconcave); wadd_(Wconcavetot, dist); wmax_(Wconcavemax, dist); }else if (mergetype == MRGconcavecoplanar) { zinc_(Zconcavecoplanar); wadd_(Wconcavecoplanartot, dist); wmax_(Wconcavecoplanarmax, dist); }else { /* MRGcoplanar */ zinc_(Zcoplanar); wadd_(Wcoplanartot, dist); wmax_(Wcoplanarmax, dist); } } } /* merge_nonconvex */ /*--------------------------------- qh_merge_pinchedvertices( apex ) merge pinched vertices in qh.vertex_mergeset to avoid qh_forcedmerges of dupridges notes: only called by qh_all_vertexmerges hull_dim >= 3 design: make vertex neighbors if necessary for each pinched vertex determine the ridges for the pinched vertex (make ridges as needed) merge the pinched vertex into the horizon vertex merge the degenerate and redundant facets that result check and resolve new dupridges */ void qh_merge_pinchedvertices(int apexpointid /* qh.newfacet_list */) { mergeT *merge, *mergeA, **mergeAp; vertexT *vertex, *vertex2; realT dist; boolT firstmerge= True; qh_vertexneighbors(); if (qh visible_list || qh newfacet_list || qh newvertex_list) { qh_fprintf(qh ferr, 6402, "qhull internal error (qh_merge_pinchedvertices): qh.visible_list (f%d), newfacet_list (f%d), or newvertex_list (v%d) not empty\n", getid_(qh visible_list), getid_(qh newfacet_list), getid_(qh newvertex_list)); qh_errexit(qh_ERRqhull, NULL, NULL); } qh visible_list= qh newfacet_list= qh facet_tail; qh newvertex_list= qh vertex_tail; qh isRenameVertex= True; /* disable duplicate ridge vertices check in qh_checkfacet */ while ((merge= qh_next_vertexmerge(/* qh.vertex_mergeset */))) { /* only one at a time from qh_getpinchedmerges */ if (qh TRACEmerge-1 == zzval_(Ztotmerge)) qhmem.IStracing= qh IStracing= qh TRACElevel; if (merge->mergetype == MRGsubridge) { zzinc_(Zpinchedvertex); trace1((qh ferr, 1050, "qh_merge_pinchedvertices: merge one of %d pinched vertices before adding apex p%d. Try to resolve duplicate ridges in newfacets\n", qh_setsize(qh vertex_mergeset)+1, apexpointid)); qh_remove_mergetype(qh vertex_mergeset, MRGsubridge); }else { zzinc_(Zpinchduplicate); if (firstmerge) trace1((qh ferr, 1056, "qh_merge_pinchedvertices: merge %d pinched vertices from dupridges in merged facets, apex p%d\n", qh_setsize(qh vertex_mergeset)+1, apexpointid)); firstmerge= False; } vertex= merge->vertex1; vertex2= merge->vertex2; dist= merge->distance; qh_memfree(merge, (int)sizeof(mergeT)); /* merge is invalidated */ qh_rename_adjacentvertex(vertex, vertex2, dist); #ifndef qh_NOtrace if (qh IStracing >= 2) { FOREACHmergeA_(qh degen_mergeset) { if (mergeA->mergetype== MRGdegen) { qh_fprintf(qh ferr, 2072, "qh_merge_pinchedvertices: merge degenerate f%d into an adjacent facet\n", mergeA->facet1->id); }else { qh_fprintf(qh ferr, 2084, "qh_merge_pinchedvertices: merge f%d into f%d mergeType %d\n", mergeA->facet1->id, mergeA->facet2->id, mergeA->mergetype); } } } #endif qh_merge_degenredundant(); /* simplicial facets with both old and new vertices */ } qh isRenameVertex= False; }/* merge_pinchedvertices */ /*--------------------------------- qh_merge_twisted( facet1, facet2 ) remove twisted ridge between facet1 into facet2 or report error returns: merges one of the facets into the best neighbor notes: a twisted ridge has opposite vertices that are convex and concave design: find best neighbors for both facets error if wide merge merge the nearest facet into its best neighbor update statistics */ void qh_merge_twisted(facetT *facet1, facetT *facet2) { facetT *neighbor2, *neighbor, *merging, *merged; vertexT *bestvertex, *bestpinched; realT dist, dist2, mindist, mindist2, maxdist, maxdist2, mintwisted, bestdist; if (qh TRACEmerge-1 == zzval_(Ztotmerge)) qhmem.IStracing= qh IStracing= qh TRACElevel; trace3((qh ferr, 3050, "qh_merge_twisted: merge #%d for twisted f%d and f%d\n", zzval_(Ztotmerge) + 1, facet1->id, facet2->id)); /* twisted */ neighbor= qh_findbestneighbor(facet1, &dist, &mindist, &maxdist); neighbor2= qh_findbestneighbor(facet2, &dist2, &mindist2, &maxdist2); mintwisted= qh_RATIOtwisted * qh ONEmerge; maximize_(mintwisted, facet1->maxoutside); maximize_(mintwisted, facet2->maxoutside); if (dist > mintwisted && dist2 > mintwisted) { bestdist= qh_vertex_bestdist2(facet1->vertices, &bestvertex, &bestpinched); if (bestdist > mintwisted) { qh_fprintf(qh ferr, 6417, "qhull precision error (qh_merge_twisted): twisted facet f%d does not contain pinched vertices. Too wide to merge into neighbor. mindist %2.2g maxdist %2.2g vertexdist %2.2g maxpinched %2.2g neighbor f%d mindist %2.2g maxdist %2.2g\n", facet1->id, mindist, maxdist, bestdist, mintwisted, facet2->id, mindist2, maxdist2); }else { qh_fprintf(qh ferr, 6418, "qhull precision error (qh_merge_twisted): twisted facet f%d with pinched vertices. Could merge vertices, but too wide to merge into neighbor. mindist %2.2g maxdist %2.2g vertexdist %2.2g neighbor f%d mindist %2.2g maxdist %2.2g\n", facet1->id, mindist, maxdist, bestdist, facet2->id, mindist2, maxdist2); } qh_errexit2(qh_ERRwide, facet1, facet2); } if (dist < dist2) { merging= facet1; merged= neighbor; }else { /* ignores qh.AVOIDold ('Q4') */ merging= facet2; merged= neighbor2; dist= dist2; mindist= mindist2; maxdist= maxdist2; } qh_mergefacet(merging, merged, MRGtwisted, &mindist, &maxdist, !qh_MERGEapex); /* caller merges qh_degenredundant */ zinc_(Ztwisted); wadd_(Wtwistedtot, dist); wmax_(Wtwistedmax, dist); } /* merge_twisted */ /*--------------------------------- qh_mergecycle( samecycle, newfacet ) merge a cycle of facets starting at samecycle into a newfacet newfacet is a horizon facet with ->normal samecycle facets are simplicial from an apex returns: initializes vertex neighbors on first merge samecycle deleted (placed on qh.visible_list) newfacet at end of qh.facet_list deleted vertices on qh.del_vertices notes: only called by qh_mergecycle_all for multiple, same cycle facets see qh_mergefacet design: make vertex neighbors if necessary make ridges for newfacet merge neighbor sets of samecycle into newfacet merge ridges of samecycle into newfacet merge vertex neighbors of samecycle into newfacet make apex of samecycle the apex of newfacet if newfacet wasn't a new facet add its vertices to qh.newvertex_list delete samecycle facets a make newfacet a newfacet */ void qh_mergecycle(facetT *samecycle, facetT *newfacet) { int traceonce= False, tracerestore= 0; vertexT *apex; #ifndef qh_NOtrace facetT *same; #endif zzinc_(Ztotmerge); if (qh REPORTfreq2 && qh POSTmerging) { if (zzval_(Ztotmerge) > qh mergereport + qh REPORTfreq2) qh_tracemerging(); } #ifndef qh_NOtrace if (qh TRACEmerge == zzval_(Ztotmerge)) qhmem.IStracing= qh IStracing= qh TRACElevel; trace2((qh ferr, 2030, "qh_mergecycle: merge #%d for facets from cycle f%d into coplanar horizon f%d\n", zzval_(Ztotmerge), samecycle->id, newfacet->id)); if (newfacet == qh tracefacet) { tracerestore= qh IStracing; qh IStracing= 4; qh_fprintf(qh ferr, 8068, "qh_mergecycle: ========= trace merge %d of samecycle %d into trace f%d, furthest is p%d\n", zzval_(Ztotmerge), samecycle->id, newfacet->id, qh furthest_id); traceonce= True; } if (qh IStracing >=4) { qh_fprintf(qh ferr, 8069, " same cycle:"); FORALLsame_cycle_(samecycle) qh_fprintf(qh ferr, 8070, " f%d", same->id); qh_fprintf(qh ferr, 8071, "\n"); } if (qh IStracing >=4) qh_errprint("MERGING CYCLE", samecycle, newfacet, NULL, NULL); #endif /* !qh_NOtrace */ if (newfacet->tricoplanar) { if (!qh TRInormals) { qh_fprintf(qh ferr, 6224, "qhull internal error (qh_mergecycle): does not work for tricoplanar facets. Use option 'Q11'\n"); qh_errexit(qh_ERRqhull, newfacet, NULL); } newfacet->tricoplanar= False; newfacet->keepcentrum= False; } if (qh CHECKfrequently) qh_checkdelridge(); if (!qh VERTEXneighbors) qh_vertexneighbors(); apex= SETfirstt_(samecycle->vertices, vertexT); qh_makeridges(newfacet); qh_mergecycle_neighbors(samecycle, newfacet); qh_mergecycle_ridges(samecycle, newfacet); qh_mergecycle_vneighbors(samecycle, newfacet); if (SETfirstt_(newfacet->vertices, vertexT) != apex) qh_setaddnth(&newfacet->vertices, 0, apex); /* apex has last id */ if (!newfacet->newfacet) qh_newvertices(newfacet->vertices); qh_mergecycle_facets(samecycle, newfacet); qh_tracemerge(samecycle, newfacet, MRGcoplanarhorizon); /* check for degen_redundant_neighbors after qh_forcedmerges() */ if (traceonce) { qh_fprintf(qh ferr, 8072, "qh_mergecycle: end of trace facet\n"); qh IStracing= tracerestore; } } /* mergecycle */ /*--------------------------------- qh_mergecycle_all( facetlist, wasmerge ) merge all samecycles of coplanar facets into horizon don't merge facets with ->mergeridge (these already have ->normal) all facets are simplicial from apex all facet->cycledone == False returns: all newfacets merged into coplanar horizon facets deleted vertices on qh.del_vertices sets wasmerge if any merge notes: called by qh_premerge calls qh_mergecycle for multiple, same cycle facets design: for each facet on facetlist skip facets with dupridges and normals check that facet is in a samecycle (->mergehorizon) if facet only member of samecycle sets vertex->delridge for all vertices except apex merge facet into horizon else mark all facets in samecycle remove facets with dupridges from samecycle merge samecycle into horizon (deletes facets from facetlist) */ void qh_mergecycle_all(facetT *facetlist, boolT *wasmerge) { facetT *facet, *same, *prev, *horizon, *newfacet; facetT *samecycle= NULL, *nextfacet, *nextsame; vertexT *apex, *vertex, **vertexp; int cycles=0, total=0, facets, nummerge, numdegen= 0; trace2((qh ferr, 2031, "qh_mergecycle_all: merge new facets into coplanar horizon facets. Bulk merge a cycle of facets with the same horizon facet\n")); for (facet=facetlist; facet && (nextfacet= facet->next); facet= nextfacet) { if (facet->normal) continue; if (!facet->mergehorizon) { qh_fprintf(qh ferr, 6225, "qhull internal error (qh_mergecycle_all): f%d without normal\n", facet->id); qh_errexit(qh_ERRqhull, facet, NULL); } horizon= SETfirstt_(facet->neighbors, facetT); if (facet->f.samecycle == facet) { if (qh TRACEmerge-1 == zzval_(Ztotmerge)) qhmem.IStracing= qh IStracing= qh TRACElevel; zinc_(Zonehorizon); /* merge distance done in qh_findhorizon */ apex= SETfirstt_(facet->vertices, vertexT); FOREACHvertex_(facet->vertices) { if (vertex != apex) vertex->delridge= True; } horizon->f.newcycle= NULL; qh_mergefacet(facet, horizon, MRGcoplanarhorizon, NULL, NULL, qh_MERGEapex); }else { samecycle= facet; facets= 0; prev= facet; for (same= facet->f.samecycle; same; /* FORALLsame_cycle_(facet) */ same= (same == facet ? NULL :nextsame)) { /* ends at facet */ nextsame= same->f.samecycle; if (same->cycledone || same->visible) qh_infiniteloop(same); same->cycledone= True; if (same->normal) { prev->f.samecycle= same->f.samecycle; /* unlink ->mergeridge */ same->f.samecycle= NULL; }else { prev= same; facets++; } } while (nextfacet && nextfacet->cycledone) /* will delete samecycle */ nextfacet= nextfacet->next; horizon->f.newcycle= NULL; qh_mergecycle(samecycle, horizon); nummerge= horizon->nummerge + facets; if (nummerge > qh_MAXnummerge) horizon->nummerge= qh_MAXnummerge; else horizon->nummerge= (short unsigned int)nummerge; /* limited to 9 bits by qh_MAXnummerge, -Wconversion */ zzinc_(Zcyclehorizon); total += facets; zzadd_(Zcyclefacettot, facets); zmax_(Zcyclefacetmax, facets); } cycles++; } if (cycles) { FORALLnew_facets { /* qh_maybe_duplicateridges postponed since qh_mergecycle_ridges deletes ridges without calling qh_delridge_merge */ if (newfacet->coplanarhorizon) { qh_test_redundant_neighbors(newfacet); qh_maybe_duplicateridges(newfacet); newfacet->coplanarhorizon= False; } } numdegen += qh_merge_degenredundant(); *wasmerge= True; trace1((qh ferr, 1013, "qh_mergecycle_all: merged %d same cycles or facets into coplanar horizons and %d degenredundant facets\n", cycles, numdegen)); } } /* mergecycle_all */ /*--------------------------------- qh_mergecycle_facets( samecycle, newfacet ) finish merge of samecycle into newfacet returns: samecycle prepended to visible_list for later deletion and partitioning each facet->f.replace == newfacet newfacet moved to end of qh.facet_list makes newfacet a newfacet (get's facet1->id if it was old) sets newfacet->newmerge clears newfacet->center (unless merging into a large facet) clears newfacet->tested and ridge->tested for facet1 adds neighboring facets to facet_mergeset if redundant or degenerate design: make newfacet a new facet and set its flags move samecycle facets to qh.visible_list for later deletion unless newfacet is large remove its centrum */ void qh_mergecycle_facets(facetT *samecycle, facetT *newfacet) { facetT *same, *next; trace4((qh ferr, 4030, "qh_mergecycle_facets: make newfacet new and samecycle deleted\n")); qh_removefacet(newfacet); /* append as a newfacet to end of qh facet_list */ qh_appendfacet(newfacet); newfacet->newfacet= True; newfacet->simplicial= False; newfacet->newmerge= True; for (same= samecycle->f.samecycle; same; same= (same == samecycle ? NULL : next)) { next= same->f.samecycle; /* reused by willdelete */ qh_willdelete(same, newfacet); } if (newfacet->center && qh_setsize(newfacet->vertices) <= qh hull_dim + qh_MAXnewcentrum) { qh_memfree(newfacet->center, qh normal_size); newfacet->center= NULL; } trace3((qh ferr, 3004, "qh_mergecycle_facets: merged facets from cycle f%d into f%d\n", samecycle->id, newfacet->id)); } /* mergecycle_facets */ /*--------------------------------- qh_mergecycle_neighbors( samecycle, newfacet ) add neighbors for samecycle facets to newfacet returns: newfacet with updated neighbors and vice-versa newfacet has ridges all neighbors of newfacet marked with qh.visit_id samecycle facets marked with qh.visit_id-1 ridges updated for simplicial neighbors of samecycle with a ridge notes: assumes newfacet not in samecycle usually, samecycle facets are new, simplicial facets without internal ridges not so if horizon facet is coplanar to two different samecycles see: qh_mergeneighbors() design: check samecycle delete neighbors from newfacet that are also in samecycle for each neighbor of a facet in samecycle if neighbor is simplicial if first visit move the neighbor relation to newfacet update facet links for its ridges else make ridges for neighbor remove samecycle reference else update neighbor sets */ void qh_mergecycle_neighbors(facetT *samecycle, facetT *newfacet) { facetT *same, *neighbor, **neighborp; int delneighbors= 0, newneighbors= 0; unsigned int samevisitid; ridgeT *ridge, **ridgep; samevisitid= ++qh visit_id; FORALLsame_cycle_(samecycle) { if (same->visitid == samevisitid || same->visible) qh_infiniteloop(samecycle); same->visitid= samevisitid; } newfacet->visitid= ++qh visit_id; trace4((qh ferr, 4031, "qh_mergecycle_neighbors: delete shared neighbors from newfacet\n")); FOREACHneighbor_(newfacet) { if (neighbor->visitid == samevisitid) { SETref_(neighbor)= NULL; /* samecycle neighbors deleted */ delneighbors++; }else neighbor->visitid= qh visit_id; } qh_setcompact(newfacet->neighbors); trace4((qh ferr, 4032, "qh_mergecycle_neighbors: update neighbors\n")); FORALLsame_cycle_(samecycle) { FOREACHneighbor_(same) { if (neighbor->visitid == samevisitid) continue; if (neighbor->simplicial) { if (neighbor->visitid != qh visit_id) { qh_setappend(&newfacet->neighbors, neighbor); qh_setreplace(neighbor->neighbors, same, newfacet); newneighbors++; neighbor->visitid= qh visit_id; FOREACHridge_(neighbor->ridges) { /* update ridge in case of qh_makeridges */ if (ridge->top == same) { ridge->top= newfacet; break; }else if (ridge->bottom == same) { ridge->bottom= newfacet; break; } } }else { qh_makeridges(neighbor); qh_setdel(neighbor->neighbors, same); /* same can't be horizon facet for neighbor */ } }else { /* non-simplicial neighbor */ qh_setdel(neighbor->neighbors, same); if (neighbor->visitid != qh visit_id) { qh_setappend(&neighbor->neighbors, newfacet); qh_setappend(&newfacet->neighbors, neighbor); neighbor->visitid= qh visit_id; newneighbors++; } } } } trace2((qh ferr, 2032, "qh_mergecycle_neighbors: deleted %d neighbors and added %d\n", delneighbors, newneighbors)); } /* mergecycle_neighbors */ /*--------------------------------- qh_mergecycle_ridges( samecycle, newfacet ) add ridges/neighbors for facets in samecycle to newfacet all new/old neighbors of newfacet marked with qh.visit_id facets in samecycle marked with qh.visit_id-1 newfacet marked with qh.visit_id returns: newfacet has merged ridges notes: ridge already updated for simplicial neighbors of samecycle with a ridge qh_checkdelridge called by qh_mergecycle see: qh_mergeridges() qh_makeridges() design: remove ridges between newfacet and samecycle for each facet in samecycle for each ridge in facet update facet pointers in ridge skip ridges processed in qh_mergecycle_neighors free ridges between newfacet and samecycle free ridges between facets of samecycle (on 2nd visit) append remaining ridges to newfacet if simplicial facet for each neighbor of facet if simplicial facet and not samecycle facet or newfacet make ridge between neighbor and newfacet */ void qh_mergecycle_ridges(facetT *samecycle, facetT *newfacet) { facetT *same, *neighbor= NULL; int numold=0, numnew=0; int neighbor_i, neighbor_n; unsigned int samevisitid; ridgeT *ridge, **ridgep; boolT toporient; void **freelistp; /* used if !qh_NOmem by qh_memfree_() */ trace4((qh ferr, 4033, "qh_mergecycle_ridges: delete shared ridges from newfacet\n")); samevisitid= qh visit_id -1; FOREACHridge_(newfacet->ridges) { neighbor= otherfacet_(ridge, newfacet); if (neighbor->visitid == samevisitid) SETref_(ridge)= NULL; /* ridge free'd below */ } qh_setcompact(newfacet->ridges); trace4((qh ferr, 4034, "qh_mergecycle_ridges: add ridges to newfacet\n")); FORALLsame_cycle_(samecycle) { FOREACHridge_(same->ridges) { if (ridge->top == same) { ridge->top= newfacet; neighbor= ridge->bottom; }else if (ridge->bottom == same) { ridge->bottom= newfacet; neighbor= ridge->top; }else if (ridge->top == newfacet || ridge->bottom == newfacet) { qh_setappend(&newfacet->ridges, ridge); numold++; /* already set by qh_mergecycle_neighbors */ continue; }else { qh_fprintf(qh ferr, 6098, "qhull internal error (qh_mergecycle_ridges): bad ridge r%d\n", ridge->id); qh_errexit(qh_ERRqhull, NULL, ridge); } if (neighbor == newfacet) { if (qh traceridge == ridge) qh traceridge= NULL; qh_setfree(&(ridge->vertices)); qh_memfree_(ridge, (int)sizeof(ridgeT), freelistp); numold++; }else if (neighbor->visitid == samevisitid) { qh_setdel(neighbor->ridges, ridge); if (qh traceridge == ridge) qh traceridge= NULL; qh_setfree(&(ridge->vertices)); qh_memfree_(ridge, (int)sizeof(ridgeT), freelistp); numold++; }else { qh_setappend(&newfacet->ridges, ridge); numold++; } } if (same->ridges) qh_settruncate(same->ridges, 0); if (!same->simplicial) continue; FOREACHneighbor_i_(same) { /* note: !newfact->simplicial */ if (neighbor->visitid != samevisitid && neighbor->simplicial) { ridge= qh_newridge(); ridge->vertices= qh_setnew_delnthsorted(same->vertices, qh hull_dim, neighbor_i, 0); toporient= (boolT)(same->toporient ^ (neighbor_i & 0x1)); if (toporient) { ridge->top= newfacet; ridge->bottom= neighbor; ridge->simplicialbot= True; }else { ridge->top= neighbor; ridge->bottom= newfacet; ridge->simplicialtop= True; } qh_setappend(&(newfacet->ridges), ridge); qh_setappend(&(neighbor->ridges), ridge); if (qh ridge_id == qh traceridge_id) qh traceridge= ridge; numnew++; } } } trace2((qh ferr, 2033, "qh_mergecycle_ridges: found %d old ridges and %d new ones\n", numold, numnew)); } /* mergecycle_ridges */ /*--------------------------------- qh_mergecycle_vneighbors( samecycle, newfacet ) create vertex neighbors for newfacet from vertices of facets in samecycle samecycle marked with visitid == qh.visit_id - 1 returns: newfacet vertices with updated neighbors marks newfacet with qh.visit_id-1 deletes vertices that are merged away sets delridge on all vertices (faster here than in mergecycle_ridges) see: qh_mergevertex_neighbors() design: for each vertex of samecycle facet set vertex->delridge delete samecycle facets from vertex neighbors append newfacet to vertex neighbors if vertex only in newfacet delete it from newfacet add it to qh.del_vertices for later deletion */ void qh_mergecycle_vneighbors(facetT *samecycle, facetT *newfacet) { facetT *neighbor, **neighborp; unsigned int mergeid; vertexT *vertex, **vertexp, *apex; setT *vertices; trace4((qh ferr, 4035, "qh_mergecycle_vneighbors: update vertex neighbors for newfacet\n")); mergeid= qh visit_id - 1; newfacet->visitid= mergeid; vertices= qh_basevertices(samecycle); /* temp */ apex= SETfirstt_(samecycle->vertices, vertexT); qh_setappend(&vertices, apex); FOREACHvertex_(vertices) { vertex->delridge= True; FOREACHneighbor_(vertex) { if (neighbor->visitid == mergeid) SETref_(neighbor)= NULL; } qh_setcompact(vertex->neighbors); qh_setappend(&vertex->neighbors, newfacet); if (!SETsecond_(vertex->neighbors)) { zinc_(Zcyclevertex); trace2((qh ferr, 2034, "qh_mergecycle_vneighbors: deleted v%d when merging cycle f%d into f%d\n", vertex->id, samecycle->id, newfacet->id)); qh_setdelsorted(newfacet->vertices, vertex); vertex->deleted= True; qh_setappend(&qh del_vertices, vertex); } } qh_settempfree(&vertices); trace3((qh ferr, 3005, "qh_mergecycle_vneighbors: merged vertices from cycle f%d into f%d\n", samecycle->id, newfacet->id)); } /* mergecycle_vneighbors */ /*--------------------------------- qh_mergefacet( facet1, facet2, mergetype, mindist, maxdist, mergeapex ) merges facet1 into facet2 mergeapex==qh_MERGEapex if merging new facet into coplanar horizon (optimizes qh_mergesimplex) returns: qh.max_outside and qh.min_vertex updated initializes vertex neighbors on first merge note: mergetype only used for logging and error reporting returns: facet2 contains facet1's vertices, neighbors, and ridges facet2 moved to end of qh.facet_list makes facet2 a newfacet sets facet2->newmerge set clears facet2->center (unless merging into a large facet) clears facet2->tested and ridge->tested for facet1 facet1 prepended to visible_list for later deletion and partitioning facet1->f.replace == facet2 adds neighboring facets to facet_mergeset if redundant or degenerate notes: when done, tests facet1 and facet2 for degenerate or redundant neighbors and dupridges mindist/maxdist may be NULL (only if both NULL) traces merge if fmax_(maxdist,-mindist) > TRACEdist see: qh_mergecycle() design: trace merge and check for degenerate simplex make ridges for both facets update qh.max_outside, qh.max_vertex, qh.min_vertex update facet2->maxoutside and keepcentrum update facet2->nummerge update tested flags for facet2 if facet1 is simplicial merge facet1 into facet2 else merge facet1's neighbors into facet2 merge facet1's ridges into facet2 merge facet1's vertices into facet2 merge facet1's vertex neighbors into facet2 add facet2's vertices to qh.new_vertexlist move facet2 to end of qh.newfacet_list unless MRGcoplanarhorizon test facet2 for redundant neighbors test facet1 for degenerate neighbors test for redundant facet2 maybe test for duplicate ridges ('Q15') move facet1 to qh.visible_list for later deletion */ void qh_mergefacet(facetT *facet1, facetT *facet2, mergeType mergetype, realT *mindist, realT *maxdist, boolT mergeapex) { boolT traceonce= False; vertexT *vertex, **vertexp; realT mintwisted, vertexdist; realT onemerge; int tracerestore=0, nummerge; const char *mergename; if(mergetype > 0 && mergetype < sizeof(mergetypes)/sizeof(char *)) mergename= mergetypes[mergetype]; else mergename= mergetypes[MRGnone]; if (facet1->tricoplanar || facet2->tricoplanar) { if (!qh TRInormals) { qh_fprintf(qh ferr, 6226, "qhull internal error (qh_mergefacet): merge f%d into f%d for mergetype %d (%s) does not work for tricoplanar facets. Use option 'Q11'\n", facet1->id, facet2->id, mergetype, mergename); qh_errexit2(qh_ERRqhull, facet1, facet2); } if (facet2->tricoplanar) { facet2->tricoplanar= False; facet2->keepcentrum= False; } } zzinc_(Ztotmerge); if (qh REPORTfreq2 && qh POSTmerging) { if (zzval_(Ztotmerge) > qh mergereport + qh REPORTfreq2) qh_tracemerging(); } #ifndef qh_NOtrace if (qh build_cnt >= qh RERUN) { if (mindist && (-*mindist > qh TRACEdist || *maxdist > qh TRACEdist)) { tracerestore= 0; qh IStracing= qh TRACElevel; traceonce= True; qh_fprintf(qh ferr, 8075, "qh_mergefacet: ========= trace wide merge #%d(%2.2g) for f%d into f%d for mergetype %d (%s), last point was p%d\n", zzval_(Ztotmerge), fmax_(-*mindist, *maxdist), facet1->id, facet2->id, mergetype, mergename, qh furthest_id); }else if (facet1 == qh tracefacet || facet2 == qh tracefacet) { tracerestore= qh IStracing; qh IStracing= 4; traceonce= True; qh_fprintf(qh ferr, 8076, "qh_mergefacet: ========= trace merge #%d for f%d into f%d for mergetype %d (%s), furthest is p%d\n", zzval_(Ztotmerge), facet1->id, facet2->id, mergetype, mergename, qh furthest_id); } } if (qh IStracing >= 2) { realT mergemin= -2; realT mergemax= -2; if (mindist) { mergemin= *mindist; mergemax= *maxdist; } qh_fprintf(qh ferr, 2081, "qh_mergefacet: #%d merge f%d into f%d for merge for mergetype %d (%s), mindist= %2.2g, maxdist= %2.2g, max_outside %2.2g\n", zzval_(Ztotmerge), facet1->id, facet2->id, mergetype, mergename, mergemin, mergemax, qh max_outside); } #endif /* !qh_NOtrace */ if(!qh ALLOWwide && mindist) { mintwisted= qh_WIDEmaxoutside * qh ONEmerge; /* same as qh_merge_twisted and qh_check_maxout (poly2) */ maximize_(mintwisted, facet1->maxoutside); maximize_(mintwisted, facet2->maxoutside); if (*maxdist > mintwisted || -*mindist > mintwisted) { vertexdist= qh_vertex_bestdist(facet1->vertices); onemerge= qh ONEmerge + qh DISTround; if (vertexdist > mintwisted) { qh_fprintf(qh ferr, 6347, "qhull precision error (qh_mergefacet): wide merge for facet f%d into f%d for mergetype %d (%s). maxdist %2.2g (%.1fx) mindist %2.2g (%.1fx) vertexdist %2.2g Allow with 'Q12' (allow-wide)\n", facet1->id, facet2->id, mergetype, mergename, *maxdist, *maxdist/onemerge, *mindist, -*mindist/onemerge, vertexdist); }else { qh_fprintf(qh ferr, 6348, "qhull precision error (qh_mergefacet): wide merge for pinched facet f%d into f%d for mergetype %d (%s). maxdist %2.2g (%.fx) mindist %2.2g (%.1fx) vertexdist %2.2g Allow with 'Q12' (allow-wide)\n", facet1->id, facet2->id, mergetype, mergename, *maxdist, *maxdist/onemerge, *mindist, -*mindist/onemerge, vertexdist); } qh_errexit2(qh_ERRwide, facet1, facet2); } } if (facet1 == facet2 || facet1->visible || facet2->visible) { qh_fprintf(qh ferr, 6099, "qhull internal error (qh_mergefacet): either f%d and f%d are the same or one is a visible facet, mergetype %d (%s)\n", facet1->id, facet2->id, mergetype, mergename); qh_errexit2(qh_ERRqhull, facet1, facet2); } if (qh num_facets - qh num_visible <= qh hull_dim + 1) { qh_fprintf(qh ferr, 6227, "qhull topology error: Only %d facets remain. The input is too degenerate or the convexity constraints are too strong.\n", qh hull_dim+1); if (qh hull_dim >= 5 && !qh MERGEexact) qh_fprintf(qh ferr, 8079, " Option 'Qx' may avoid this problem.\n"); qh_errexit(qh_ERRtopology, NULL, NULL); } if (!qh VERTEXneighbors) qh_vertexneighbors(); qh_makeridges(facet1); qh_makeridges(facet2); if (qh IStracing >=4) qh_errprint("MERGING", facet1, facet2, NULL, NULL); if (mindist) { maximize_(qh max_outside, *maxdist); maximize_(qh max_vertex, *maxdist); #if qh_MAXoutside maximize_(facet2->maxoutside, *maxdist); #endif minimize_(qh min_vertex, *mindist); if (!facet2->keepcentrum && (*maxdist > qh WIDEfacet || *mindist < -qh WIDEfacet)) { facet2->keepcentrum= True; zinc_(Zwidefacet); } } nummerge= facet1->nummerge + facet2->nummerge + 1; if (nummerge >= qh_MAXnummerge) facet2->nummerge= qh_MAXnummerge; else facet2->nummerge= (short unsigned int)nummerge; /* limited to 9 bits by qh_MAXnummerge, -Wconversion */ facet2->newmerge= True; facet2->dupridge= False; qh_updatetested(facet1, facet2); if (qh hull_dim > 2 && qh_setsize(facet1->vertices) == qh hull_dim) qh_mergesimplex(facet1, facet2, mergeapex); else { qh vertex_visit++; FOREACHvertex_(facet2->vertices) vertex->visitid= qh vertex_visit; if (qh hull_dim == 2) qh_mergefacet2d(facet1, facet2); else { qh_mergeneighbors(facet1, facet2); qh_mergevertices(facet1->vertices, &facet2->vertices); } qh_mergeridges(facet1, facet2); qh_mergevertex_neighbors(facet1, facet2); if (!facet2->newfacet) qh_newvertices(facet2->vertices); } if (facet2->coplanarhorizon) { zinc_(Zmergeintocoplanar); }else if (!facet2->newfacet) { zinc_(Zmergeintohorizon); }else if (!facet1->newfacet && facet2->newfacet) { zinc_(Zmergehorizon); }else { zinc_(Zmergenew); } qh_removefacet(facet2); /* append as a newfacet to end of qh facet_list */ qh_appendfacet(facet2); facet2->newfacet= True; facet2->tested= False; qh_tracemerge(facet1, facet2, mergetype); if (traceonce) { qh_fprintf(qh ferr, 8080, "qh_mergefacet: end of wide tracing\n"); qh IStracing= tracerestore; } if (mergetype != MRGcoplanarhorizon) { trace3((qh ferr, 3076, "qh_mergefacet: check f%d and f%d for redundant and degenerate neighbors\n", facet1->id, facet2->id)); qh_test_redundant_neighbors(facet2); qh_test_degen_neighbors(facet1); /* after qh_test_redundant_neighbors since MRGdegen more difficult than MRGredundant and before qh_willdelete which clears facet1.neighbors */ qh_degen_redundant_facet(facet2); /* may occur in qh_merge_pinchedvertices, e.g., rbox 175 C3,2e-13 D4 t1545228104 | qhull d */ qh_maybe_duplicateridges(facet2); } qh_willdelete(facet1, facet2); } /* mergefacet */ /*--------------------------------- qh_mergefacet2d( facet1, facet2 ) in 2d, merges neighbors and vertices of facet1 into facet2 returns: build ridges for neighbors if necessary facet2 looks like a simplicial facet except for centrum, ridges neighbors are opposite the corresponding vertex maintains orientation of facet2 notes: qh_mergefacet() retains non-simplicial structures they are not needed in 2d, but later routines may use them preserves qh.vertex_visit for qh_mergevertex_neighbors() design: get vertices and neighbors determine new vertices and neighbors set new vertices and neighbors and adjust orientation make ridges for new neighbor if needed */ void qh_mergefacet2d(facetT *facet1, facetT *facet2) { vertexT *vertex1A, *vertex1B, *vertex2A, *vertex2B, *vertexA, *vertexB; facetT *neighbor1A, *neighbor1B, *neighbor2A, *neighbor2B, *neighborA, *neighborB; vertex1A= SETfirstt_(facet1->vertices, vertexT); vertex1B= SETsecondt_(facet1->vertices, vertexT); vertex2A= SETfirstt_(facet2->vertices, vertexT); vertex2B= SETsecondt_(facet2->vertices, vertexT); neighbor1A= SETfirstt_(facet1->neighbors, facetT); neighbor1B= SETsecondt_(facet1->neighbors, facetT); neighbor2A= SETfirstt_(facet2->neighbors, facetT); neighbor2B= SETsecondt_(facet2->neighbors, facetT); if (vertex1A == vertex2A) { vertexA= vertex1B; vertexB= vertex2B; neighborA= neighbor2A; neighborB= neighbor1A; }else if (vertex1A == vertex2B) { vertexA= vertex1B; vertexB= vertex2A; neighborA= neighbor2B; neighborB= neighbor1A; }else if (vertex1B == vertex2A) { vertexA= vertex1A; vertexB= vertex2B; neighborA= neighbor2A; neighborB= neighbor1B; }else { /* 1B == 2B */ vertexA= vertex1A; vertexB= vertex2A; neighborA= neighbor2B; neighborB= neighbor1B; } /* vertexB always from facet2, neighborB always from facet1 */ if (vertexA->id > vertexB->id) { SETfirst_(facet2->vertices)= vertexA; SETsecond_(facet2->vertices)= vertexB; if (vertexB == vertex2A) facet2->toporient= !facet2->toporient; SETfirst_(facet2->neighbors)= neighborA; SETsecond_(facet2->neighbors)= neighborB; }else { SETfirst_(facet2->vertices)= vertexB; SETsecond_(facet2->vertices)= vertexA; if (vertexB == vertex2B) facet2->toporient= !facet2->toporient; SETfirst_(facet2->neighbors)= neighborB; SETsecond_(facet2->neighbors)= neighborA; } /* qh_makeridges not needed since neighborB is not degenerate */ qh_setreplace(neighborB->neighbors, facet1, facet2); trace4((qh ferr, 4036, "qh_mergefacet2d: merged v%d and neighbor f%d of f%d into f%d\n", vertexA->id, neighborB->id, facet1->id, facet2->id)); } /* mergefacet2d */ /*--------------------------------- qh_mergeneighbors( facet1, facet2 ) merges the neighbors of facet1 into facet2 notes: only called by qh_mergefacet qh.hull_dim >= 3 see qh_mergecycle_neighbors design: for each neighbor of facet1 if neighbor is also a neighbor of facet2 if neighbor is simplicial make ridges for later deletion as a degenerate facet update its neighbor set else move the neighbor relation to facet2 remove the neighbor relation for facet1 and facet2 */ void qh_mergeneighbors(facetT *facet1, facetT *facet2) { facetT *neighbor, **neighborp; trace4((qh ferr, 4037, "qh_mergeneighbors: merge neighbors of f%d and f%d\n", facet1->id, facet2->id)); qh visit_id++; FOREACHneighbor_(facet2) { neighbor->visitid= qh visit_id; } FOREACHneighbor_(facet1) { if (neighbor->visitid == qh visit_id) { if (neighbor->simplicial) /* is degen, needs ridges */ qh_makeridges(neighbor); if (SETfirstt_(neighbor->neighbors, facetT) != facet1) /*keep newfacet->horizon*/ qh_setdel(neighbor->neighbors, facet1); else { qh_setdel(neighbor->neighbors, facet2); qh_setreplace(neighbor->neighbors, facet1, facet2); } }else if (neighbor != facet2) { qh_setappend(&(facet2->neighbors), neighbor); qh_setreplace(neighbor->neighbors, facet1, facet2); } } qh_setdel(facet1->neighbors, facet2); /* here for makeridges */ qh_setdel(facet2->neighbors, facet1); } /* mergeneighbors */ /*--------------------------------- qh_mergeridges( facet1, facet2 ) merges the ridge set of facet1 into facet2 returns: may delete all ridges for a vertex sets vertex->delridge on deleted ridges see: qh_mergecycle_ridges() design: delete ridges between facet1 and facet2 mark (delridge) vertices on these ridges for later testing for each remaining ridge rename facet1 to facet2 */ void qh_mergeridges(facetT *facet1, facetT *facet2) { ridgeT *ridge, **ridgep; trace4((qh ferr, 4038, "qh_mergeridges: merge ridges of f%d into f%d\n", facet1->id, facet2->id)); FOREACHridge_(facet2->ridges) { if ((ridge->top == facet1) || (ridge->bottom == facet1)) { /* ridge.nonconvex is irrelevant due to merge */ qh_delridge_merge(ridge); /* expensive in high-d, could rebuild */ ridgep--; /* deleted this ridge, repeat with next ridge*/ } } FOREACHridge_(facet1->ridges) { if (ridge->top == facet1) { ridge->top= facet2; ridge->simplicialtop= False; }else { /* ridge.bottom is facet1 */ ridge->bottom= facet2; ridge->simplicialbot= False; } qh_setappend(&(facet2->ridges), ridge); } } /* mergeridges */ /*--------------------------------- qh_mergesimplex( facet1, facet2, mergeapex ) merge simplicial facet1 into facet2 mergeapex==qh_MERGEapex if merging samecycle into horizon facet vertex id is latest (most recently created) facet1 may be contained in facet2 ridges exist for both facets returns: facet2 with updated vertices, ridges, neighbors updated neighbors for facet1's vertices facet1 not deleted sets vertex->delridge on deleted ridges notes: special case code since this is the most common merge called from qh_mergefacet() design: if qh_MERGEapex add vertices of facet2 to qh.new_vertexlist if necessary add apex to facet2 else for each ridge between facet1 and facet2 set vertex->delridge determine the apex for facet1 (i.e., vertex to be merged) unless apex already in facet2 insert apex into vertices for facet2 add vertices of facet2 to qh.new_vertexlist if necessary add apex to qh.new_vertexlist if necessary for each vertex of facet1 if apex rename facet1 to facet2 in its vertex neighbors else delete facet1 from vertex neighbors if only in facet2 add vertex to qh.del_vertices for later deletion for each ridge of facet1 delete ridges between facet1 and facet2 append other ridges to facet2 after renaming facet to facet2 */ void qh_mergesimplex(facetT *facet1, facetT *facet2, boolT mergeapex) { vertexT *vertex, **vertexp, *opposite; ridgeT *ridge, **ridgep; boolT isnew= False; facetT *neighbor, **neighborp, *otherfacet; if (mergeapex) { opposite= SETfirstt_(facet1->vertices, vertexT); /* apex is opposite facet2. It has the last vertex id */ trace4((qh ferr, 4086, "qh_mergesimplex: merge apex v%d of f%d into facet f%d\n", opposite->id, facet1->id, facet2->id)); if (!facet2->newfacet) qh_newvertices(facet2->vertices); /* apex, the first vertex, is already new */ if (SETfirstt_(facet2->vertices, vertexT) != opposite) { qh_setaddnth(&facet2->vertices, 0, opposite); isnew= True; } }else { zinc_(Zmergesimplex); FOREACHvertex_(facet1->vertices) vertex->seen= False; FOREACHridge_(facet1->ridges) { if (otherfacet_(ridge, facet1) == facet2) { FOREACHvertex_(ridge->vertices) { vertex->seen= True; vertex->delridge= True; } break; } } FOREACHvertex_(facet1->vertices) { if (!vertex->seen) break; /* must occur */ } opposite= vertex; trace4((qh ferr, 4039, "qh_mergesimplex: merge opposite v%d of f%d into facet f%d\n", opposite->id, facet1->id, facet2->id)); isnew= qh_addfacetvertex(facet2, opposite); if (!facet2->newfacet) qh_newvertices(facet2->vertices); else if (!opposite->newfacet) { qh_removevertex(opposite); qh_appendvertex(opposite); } } trace4((qh ferr, 4040, "qh_mergesimplex: update vertex neighbors of f%d\n", facet1->id)); FOREACHvertex_(facet1->vertices) { if (vertex == opposite && isnew) qh_setreplace(vertex->neighbors, facet1, facet2); else { qh_setdel(vertex->neighbors, facet1); if (!SETsecond_(vertex->neighbors)) qh_mergevertex_del(vertex, facet1, facet2); } } trace4((qh ferr, 4041, "qh_mergesimplex: merge ridges and neighbors of f%d into f%d\n", facet1->id, facet2->id)); qh visit_id++; FOREACHneighbor_(facet2) neighbor->visitid= qh visit_id; FOREACHridge_(facet1->ridges) { otherfacet= otherfacet_(ridge, facet1); if (otherfacet == facet2) { /* ridge.nonconvex is irrelevant due to merge */ qh_delridge_merge(ridge); /* expensive in high-d, could rebuild */ ridgep--; /* deleted this ridge, repeat with next ridge*/ qh_setdel(facet2->neighbors, facet1); /* a simplicial facet may have duplicate neighbors, need to delete each one */ }else if (otherfacet->dupridge && !qh_setin(otherfacet->neighbors, facet1)) { qh_fprintf(qh ferr, 6356, "qhull topology error (qh_mergesimplex): f%d is a dupridge of f%d, cannot merge f%d into f%d\n", facet1->id, otherfacet->id, facet1->id, facet2->id); qh_errexit2(qh_ERRqhull, facet1, otherfacet); }else { trace4((qh ferr, 4059, "qh_mergesimplex: move r%d with f%d to f%d, new neighbor? %d, maybe horizon? %d\n", ridge->id, otherfacet->id, facet2->id, (otherfacet->visitid != qh visit_id), (SETfirstt_(otherfacet->neighbors, facetT) == facet1))); qh_setappend(&facet2->ridges, ridge); if (otherfacet->visitid != qh visit_id) { qh_setappend(&facet2->neighbors, otherfacet); qh_setreplace(otherfacet->neighbors, facet1, facet2); otherfacet->visitid= qh visit_id; }else { if (otherfacet->simplicial) /* is degen, needs ridges */ qh_makeridges(otherfacet); if (SETfirstt_(otherfacet->neighbors, facetT) == facet1) { /* keep new, otherfacet->neighbors->horizon */ qh_setdel(otherfacet->neighbors, facet2); qh_setreplace(otherfacet->neighbors, facet1, facet2); }else { /* facet2 is already a neighbor of otherfacet, by f.visitid */ qh_setdel(otherfacet->neighbors, facet1); } } if (ridge->top == facet1) { /* wait until after qh_makeridges */ ridge->top= facet2; ridge->simplicialtop= False; }else { ridge->bottom= facet2; ridge->simplicialbot= False; } } } trace3((qh ferr, 3006, "qh_mergesimplex: merged simplex f%d v%d into facet f%d\n", facet1->id, opposite->id, facet2->id)); } /* mergesimplex */ /*--------------------------------- qh_mergevertex_del( vertex, facet1, facet2 ) delete a vertex because of merging facet1 into facet2 returns: deletes vertex from facet2 adds vertex to qh.del_vertices for later deletion */ void qh_mergevertex_del(vertexT *vertex, facetT *facet1, facetT *facet2) { zinc_(Zmergevertex); trace2((qh ferr, 2035, "qh_mergevertex_del: deleted v%d when merging f%d into f%d\n", vertex->id, facet1->id, facet2->id)); qh_setdelsorted(facet2->vertices, vertex); vertex->deleted= True; qh_setappend(&qh del_vertices, vertex); } /* mergevertex_del */ /*--------------------------------- qh_mergevertex_neighbors( facet1, facet2 ) merge the vertex neighbors of facet1 to facet2 returns: if vertex is current qh.vertex_visit deletes facet1 from vertex->neighbors else renames facet1 to facet2 in vertex->neighbors deletes vertices if only one neighbor notes: assumes vertex neighbor sets are good */ void qh_mergevertex_neighbors(facetT *facet1, facetT *facet2) { vertexT *vertex, **vertexp; trace4((qh ferr, 4042, "qh_mergevertex_neighbors: merge vertex neighborset for f%d into f%d\n", facet1->id, facet2->id)); if (qh tracevertex) { qh_fprintf(qh ferr, 8081, "qh_mergevertex_neighbors: of f%d into f%d at furthest p%d f0= %p\n", facet1->id, facet2->id, qh furthest_id, qh tracevertex->neighbors->e[0].p); qh_errprint("TRACE", NULL, NULL, NULL, qh tracevertex); } FOREACHvertex_(facet1->vertices) { if (vertex->visitid != qh vertex_visit) qh_setreplace(vertex->neighbors, facet1, facet2); else { qh_setdel(vertex->neighbors, facet1); if (!SETsecond_(vertex->neighbors)) qh_mergevertex_del(vertex, facet1, facet2); } } if (qh tracevertex) qh_errprint("TRACE", NULL, NULL, NULL, qh tracevertex); } /* mergevertex_neighbors */ /*--------------------------------- qh_mergevertices( vertices1, vertices2 ) merges the vertex set of facet1 into facet2 returns: replaces vertices2 with merged set preserves vertex_visit for qh_mergevertex_neighbors updates qh.newvertex_list design: create a merged set of both vertices (in inverse id order) */ void qh_mergevertices(setT *vertices1, setT **vertices2) { int newsize= qh_setsize(vertices1)+qh_setsize(*vertices2) - qh hull_dim + 1; setT *mergedvertices; vertexT *vertex, **vertexp, **vertex2= SETaddr_(*vertices2, vertexT); mergedvertices= qh_settemp(newsize); FOREACHvertex_(vertices1) { if (!*vertex2 || vertex->id > (*vertex2)->id) qh_setappend(&mergedvertices, vertex); else { while (*vertex2 && (*vertex2)->id > vertex->id) qh_setappend(&mergedvertices, *vertex2++); if (!*vertex2 || (*vertex2)->id < vertex->id) qh_setappend(&mergedvertices, vertex); else qh_setappend(&mergedvertices, *vertex2++); } } while (*vertex2) qh_setappend(&mergedvertices, *vertex2++); if (newsize < qh_setsize(mergedvertices)) { qh_fprintf(qh ferr, 6100, "qhull internal error (qh_mergevertices): facets did not share a ridge\n"); qh_errexit(qh_ERRqhull, NULL, NULL); } qh_setfree(vertices2); *vertices2= mergedvertices; qh_settemppop(); } /* mergevertices */ /*--------------------------------- qh_neighbor_intersections( vertex ) return intersection of all vertices in vertex->neighbors except for vertex returns: returns temporary set of vertices does not include vertex NULL if a neighbor is simplicial NULL if empty set notes: only called by qh_redundant_vertex for qh_reducevertices so f.vertices does not contain extraneous vertices that are not in f.ridges used for renaming vertices design: initialize the intersection set with vertices of the first two neighbors delete vertex from the intersection for each remaining neighbor intersect its vertex set with the intersection set return NULL if empty return the intersection set */ setT *qh_neighbor_intersections(vertexT *vertex) { facetT *neighbor, **neighborp, *neighborA, *neighborB; setT *intersect; int neighbor_i, neighbor_n; FOREACHneighbor_(vertex) { if (neighbor->simplicial) return NULL; } neighborA= SETfirstt_(vertex->neighbors, facetT); neighborB= SETsecondt_(vertex->neighbors, facetT); zinc_(Zintersectnum); if (!neighborA) return NULL; if (!neighborB) intersect= qh_setcopy(neighborA->vertices, 0); else intersect= qh_vertexintersect_new(neighborA->vertices, neighborB->vertices); qh_settemppush(intersect); qh_setdelsorted(intersect, vertex); FOREACHneighbor_i_(vertex) { if (neighbor_i >= 2) { zinc_(Zintersectnum); qh_vertexintersect(&intersect, neighbor->vertices); if (!SETfirst_(intersect)) { zinc_(Zintersectfail); qh_settempfree(&intersect); return NULL; } } } trace3((qh ferr, 3007, "qh_neighbor_intersections: %d vertices in neighbor intersection of v%d\n", qh_setsize(intersect), vertex->id)); return intersect; } /* neighbor_intersections */ /*--------------------------------- qh_neighbor_vertices( vertex ) return neighboring vertices for a vertex (not in subridge) assumes vertices have full vertex->neighbors returns: temporary set of vertices notes: updates qh.visit_id and qh.vertex_visit similar to qh_vertexridges */ setT *qh_neighbor_vertices(vertexT *vertexA, setT *subridge) { facetT *neighbor, **neighborp; vertexT *vertex, **vertexp; setT *vertices= qh_settemp(qh TEMPsize); qh visit_id++; FOREACHneighbor_(vertexA) neighbor->visitid= qh visit_id; qh vertex_visit++; vertexA->visitid= qh vertex_visit; FOREACHvertex_(subridge) { vertex->visitid= qh vertex_visit; } FOREACHneighbor_(vertexA) { if (*neighborp) /* no new ridges in last neighbor */ qh_neighbor_vertices_facet(vertexA, neighbor, &vertices); } trace3((qh ferr, 3035, "qh_neighbor_vertices: %d non-subridge, vertex neighbors for v%d\n", qh_setsize(vertices), vertexA->id)); return vertices; } /* neighbor_vertices */ /*--------------------------------- qh_neighbor_vertices_facet( vertex, facet, vertices ) add neighboring vertices on ridges for vertex in facet neighbor->visitid==qh.visit_id if it hasn't been visited v.visitid==qh.vertex_visit if it is already in vertices returns: vertices updated sets facet->visitid to qh.visit_id-1 notes: only called by qh_neighbor_vertices similar to qh_vertexridges_facet design: for each ridge of facet if ridge of visited neighbor (i.e., unprocessed) if vertex in ridge append unprocessed vertices of ridge mark facet processed */ void qh_neighbor_vertices_facet(vertexT *vertexA, facetT *facet, setT **vertices) { ridgeT *ridge, **ridgep; facetT *neighbor; vertexT *second, *last, *vertex, **vertexp; int last_i= qh hull_dim-2, count= 0; boolT isridge; if (facet->simplicial) { FOREACHvertex_(facet->vertices) { if (vertex->visitid != qh vertex_visit) { vertex->visitid= qh vertex_visit; qh_setappend(vertices, vertex); count++; } } }else { FOREACHridge_(facet->ridges) { neighbor= otherfacet_(ridge, facet); if (neighbor->visitid == qh visit_id) { isridge= False; if (SETfirst_(ridge->vertices) == vertexA) { isridge= True; }else if (last_i > 2) { second= SETsecondt_(ridge->vertices, vertexT); last= SETelemt_(ridge->vertices, last_i, vertexT); if (second->id >= vertexA->id && last->id <= vertexA->id) { /* vertices inverse sorted by id */ if (second == vertexA || last == vertexA) isridge= True; else if (qh_setin(ridge->vertices, vertexA)) isridge= True; } }else if (SETelem_(ridge->vertices, last_i) == vertexA) { isridge= True; }else if (last_i > 1 && SETsecond_(ridge->vertices) == vertexA) { isridge= True; } if (isridge) { FOREACHvertex_(ridge->vertices) { if (vertex->visitid != qh vertex_visit) { vertex->visitid= qh vertex_visit; qh_setappend(vertices, vertex); count++; } } } } } } facet->visitid= qh visit_id-1; if (count) { trace4((qh ferr, 4079, "qh_neighbor_vertices_facet: found %d vertex neighbors for v%d in f%d (simplicial? %d)\n", count, vertexA->id, facet->id, facet->simplicial)); } } /* neighbor_vertices_facet */ /*--------------------------------- qh_newvertices( vertices ) add vertices to end of qh.vertex_list (marks as new vertices) returns: vertices on qh.newvertex_list vertex->newfacet set */ void qh_newvertices(setT *vertices) { vertexT *vertex, **vertexp; FOREACHvertex_(vertices) { if (!vertex->newfacet) { qh_removevertex(vertex); qh_appendvertex(vertex); } } } /* newvertices */ /*--------------------------------- qh_next_vertexmerge( ) return next vertex merge from qh.vertex_mergeset returns: vertex merge either MRGvertices or MRGsubridge drops merges of deleted vertices notes: called from qh_merge_pinchedvertices */ mergeT *qh_next_vertexmerge(void /* qh.vertex_mergeset */) { mergeT *merge; int merge_i, merge_n, best_i= -1; realT bestdist= REALmax; FOREACHmerge_i_(qh vertex_mergeset) { if (!merge->vertex1 || !merge->vertex2) { qh_fprintf(qh ferr, 6299, "qhull internal error (qh_next_vertexmerge): expecting two vertices for vertex merge. Got v%d v%d and optional f%d\n", getid_(merge->vertex1), getid_(merge->vertex2), getid_(merge->facet1)); qh_errexit(qh_ERRqhull, merge->facet1, NULL); } if (merge->vertex1->deleted || merge->vertex2->deleted) { trace3((qh ferr, 3030, "qh_next_vertexmerge: drop merge of v%d (del? %d) into v%d (del? %d) due to deleted vertex of r%d and r%d\n", merge->vertex1->id, merge->vertex1->deleted, merge->vertex2->id, merge->vertex2->deleted, getid_(merge->ridge1), getid_(merge->ridge2))); qh_drop_mergevertex(merge); qh_setdelnth(qh vertex_mergeset, merge_i); merge_i--; merge_n--; qh_memfree(merge, (int)sizeof(mergeT)); }else if (merge->distance < bestdist) { bestdist= merge->distance; best_i= merge_i; } } merge= NULL; if (best_i >= 0) { merge= SETelemt_(qh vertex_mergeset, best_i, mergeT); if (bestdist/qh ONEmerge > qh_WIDEpinched) { if (merge->mergetype==MRGvertices) { if (merge->ridge1->top == merge->ridge2->bottom && merge->ridge1->bottom == merge->ridge2->top) qh_fprintf(qh ferr, 6391, "qhull topology error (qh_next_vertexmerge): no nearly adjacent vertices to resolve opposite oriented ridges r%d and r%d in f%d and f%d. Nearest v%d and v%d dist %2.2g (%.1fx)\n", merge->ridge1->id, merge->ridge2->id, merge->ridge1->top->id, merge->ridge1->bottom->id, merge->vertex1->id, merge->vertex2->id, bestdist, bestdist/qh ONEmerge); else qh_fprintf(qh ferr, 6381, "qhull topology error (qh_next_vertexmerge): no nearly adjacent vertices to resolve duplicate ridges r%d and r%d. Nearest v%d and v%d dist %2.2g (%.1fx)\n", merge->ridge1->id, merge->ridge2->id, merge->vertex1->id, merge->vertex2->id, bestdist, bestdist/qh ONEmerge); }else { qh_fprintf(qh ferr, 6208, "qhull topology error (qh_next_vertexmerge): no nearly adjacent vertices to resolve dupridge. Nearest v%d and v%d dist %2.2g (%.1fx)\n", merge->vertex1->id, merge->vertex2->id, bestdist, bestdist/qh ONEmerge); } /* it may be possible to find a different vertex, after other vertex merges have occurred */ qh_errexit(qh_ERRtopology, NULL, merge->ridge1); } qh_setdelnth(qh vertex_mergeset, best_i); } return merge; } /* next_vertexmerge */ /*--------------------------------- qh_opposite_horizonfacet( merge, opposite ) return horizon facet for one of the merge facets, and its opposite vertex across the ridge assumes either facet1 or facet2 of merge is 'mergehorizon' assumes both facets are simplicial facets on qh.new_facetlist returns: horizon facet and opposite vertex notes: called by qh_getpinchedmerges */ facetT *qh_opposite_horizonfacet(mergeT *merge, vertexT **opposite) { facetT *facet, *horizon, *otherfacet; int neighbor_i; if (!merge->facet1->simplicial || !merge->facet2->simplicial || (!merge->facet1->mergehorizon && !merge->facet2->mergehorizon)) { qh_fprintf(qh ferr, 6273, "qhull internal error (qh_opposite_horizonfacet): expecting merge of simplicial facets, at least one of which is mergehorizon. Either simplicial or mergehorizon is wrong\n"); qh_errexit2(qh_ERRqhull, merge->facet1, merge->facet2); } if (merge->facet1->mergehorizon) { facet= merge->facet1; otherfacet= merge->facet2; }else { facet= merge->facet2; otherfacet= merge->facet1; } horizon= SETfirstt_(facet->neighbors, facetT); neighbor_i= qh_setindex(otherfacet->neighbors, facet); if (neighbor_i==-1) neighbor_i= qh_setindex(otherfacet->neighbors, qh_MERGEridge); if (neighbor_i==-1) { qh_fprintf(qh ferr, 6238, "qhull internal error (qh_opposite_horizonfacet): merge facet f%d not connected to mergehorizon f%d\n", otherfacet->id, facet->id); qh_errexit2(qh_ERRqhull, otherfacet, facet); } *opposite= SETelemt_(otherfacet->vertices, neighbor_i, vertexT); return horizon; } /* opposite_horizonfacet */ /*--------------------------------- qh_reducevertices( ) reduce extra vertices, shared vertices, and redundant vertices facet->newmerge is set if merged since last call vertex->delridge is set if vertex was on a deleted ridge if !qh.MERGEvertices, only removes extra vertices returns: True if also merged degen_redundant facets vertices are renamed if possible clears facet->newmerge and vertex->delridge notes: called by qh_all_merges and qh_postmerge ignored if 2-d design: merge any degenerate or redundant facets repeat until no more degenerate or redundant facets for each newly merged facet remove extra vertices if qh.MERGEvertices for each newly merged facet for each vertex if vertex was on a deleted ridge rename vertex if it is shared for each new, undeleted vertex remove delridge flag if vertex is redundant merge degenerate or redundant facets */ boolT qh_reducevertices(void) { int numshare=0, numrename= 0; boolT degenredun= False; facetT *newfacet; vertexT *vertex, **vertexp; if (qh hull_dim == 2) return False; trace2((qh ferr, 2101, "qh_reducevertices: reduce extra vertices, shared vertices, and redundant vertices\n")); if (qh_merge_degenredundant()) degenredun= True; LABELrestart: FORALLnew_facets { if (newfacet->newmerge) { if (!qh MERGEvertices) newfacet->newmerge= False; if (qh_remove_extravertices(newfacet)) { qh_degen_redundant_facet(newfacet); if (qh_merge_degenredundant()) { degenredun= True; goto LABELrestart; } } } } if (!qh MERGEvertices) return False; FORALLnew_facets { if (newfacet->newmerge) { newfacet->newmerge= False; FOREACHvertex_(newfacet->vertices) { if (vertex->delridge) { if (qh_rename_sharedvertex(vertex, newfacet)) { numshare++; if (qh_merge_degenredundant()) { degenredun= True; goto LABELrestart; } vertexp--; /* repeat since deleted vertex */ } } } } } FORALLvertex_(qh newvertex_list) { if (vertex->delridge && !vertex->deleted) { vertex->delridge= False; if (qh hull_dim >= 4 && qh_redundant_vertex(vertex)) { numrename++; if (qh_merge_degenredundant()) { degenredun= True; goto LABELrestart; } } } } trace1((qh ferr, 1014, "qh_reducevertices: renamed %d shared vertices and %d redundant vertices. Degen? %d\n", numshare, numrename, degenredun)); return degenredun; } /* reducevertices */ /*--------------------------------- qh_redundant_vertex( vertex ) rename a redundant vertex if qh_find_newvertex succeeds assumes vertices have full vertex->neighbors returns: if find a replacement vertex returns new vertex qh_renamevertex sets vertex->deleted for redundant vertex notes: only called by qh_reducevertices for vertex->delridge and hull_dim >= 4 may add degenerate facets to qh.facet_mergeset doesn't change vertex->neighbors or create redundant facets design: intersect vertices of all facet neighbors of vertex determine ridges for these vertices if find a new vertex for vertex among these ridges and vertices rename vertex to the new vertex */ vertexT *qh_redundant_vertex(vertexT *vertex) { vertexT *newvertex= NULL; setT *vertices, *ridges; trace3((qh ferr, 3008, "qh_redundant_vertex: check if v%d from a deleted ridge can be renamed\n", vertex->id)); if ((vertices= qh_neighbor_intersections(vertex))) { ridges= qh_vertexridges(vertex, !qh_ALL); if ((newvertex= qh_find_newvertex(vertex, vertices, ridges))) { zinc_(Zrenameall); qh_renamevertex(vertex, newvertex, ridges, NULL, NULL); /* ridges invalidated */ } qh_settempfree(&ridges); qh_settempfree(&vertices); } return newvertex; } /* redundant_vertex */ /*--------------------------------- qh_remove_extravertices( facet ) remove extra vertices from non-simplicial facets returns: returns True if it finds them deletes facet from vertex neighbors facet may be redundant (test with qh_degen_redundant) notes: called by qh_renamevertex and qh_reducevertices a merge (qh_reducevertices) or qh_renamevertex may drop all ridges for a vertex in a facet design: for each vertex in facet if vertex not in a ridge (i.e., no longer used) delete vertex from facet delete facet from vertex's neighbors unless vertex in another facet add vertex to qh.del_vertices for later deletion */ boolT qh_remove_extravertices(facetT *facet) { ridgeT *ridge, **ridgep; vertexT *vertex, **vertexp; boolT foundrem= False; if (facet->simplicial) { return False; } trace4((qh ferr, 4043, "qh_remove_extravertices: test non-simplicial f%d for extra vertices\n", facet->id)); FOREACHvertex_(facet->vertices) vertex->seen= False; FOREACHridge_(facet->ridges) { FOREACHvertex_(ridge->vertices) vertex->seen= True; } FOREACHvertex_(facet->vertices) { if (!vertex->seen) { foundrem= True; zinc_(Zremvertex); qh_setdelsorted(facet->vertices, vertex); qh_setdel(vertex->neighbors, facet); if (!qh_setsize(vertex->neighbors)) { vertex->deleted= True; qh_setappend(&qh del_vertices, vertex); zinc_(Zremvertexdel); trace2((qh ferr, 2036, "qh_remove_extravertices: v%d deleted because it's lost all ridges\n", vertex->id)); }else trace3((qh ferr, 3009, "qh_remove_extravertices: v%d removed from f%d because it's lost all ridges\n", vertex->id, facet->id)); vertexp--; /*repeat*/ } } return foundrem; } /* remove_extravertices */ /*--------------------------------- qh_remove_mergetype( mergeset, mergetype ) Remove mergetype merges from mergeset notes: Does not preserve order */ void qh_remove_mergetype(setT *mergeset, mergeType type) { mergeT *merge; int merge_i, merge_n; FOREACHmerge_i_(mergeset) { if (merge->mergetype == type) { trace3((qh ferr, 3037, "qh_remove_mergetype: remove merge f%d f%d v%d v%d r%d r%d dist %2.2g type %d", getid_(merge->facet1), getid_(merge->facet2), getid_(merge->vertex1), getid_(merge->vertex2), getid_(merge->ridge1), getid_(merge->ridge2), merge->distance, type)); qh_setdelnth(mergeset, merge_i); merge_i--; merge_n--; /* repeat with next merge */ } } } /* remove_mergetype */ /*--------------------------------- qh_rename_adjacentvertex( oldvertex, newvertex ) renames oldvertex as newvertex. Must be adjacent (i.e., in the same subridge) no-op if either vertex is deleted notes: called from qh_merge_pinchedvertices design: for all neighbors of oldvertex if simplicial, rename oldvertex to newvertex and drop if degenerate if needed, add oldvertex neighbor to newvertex determine ridges for oldvertex rename oldvertex as newvertex in ridges (qh_renamevertex) */ void qh_rename_adjacentvertex(vertexT *oldvertex, vertexT *newvertex, realT dist) { setT *ridges; facetT *neighbor, **neighborp, *maxfacet= NULL; ridgeT *ridge, **ridgep; boolT istrace= False; int oldsize= qh_setsize(oldvertex->neighbors); int newsize= qh_setsize(newvertex->neighbors); coordT maxdist2= -REALmax, dist2; if (qh IStracing >= 4 || oldvertex->id == qh tracevertex_id || newvertex->id == qh tracevertex_id) { istrace= True; } zzinc_(Ztotmerge); if (istrace) { qh_fprintf(qh ferr, 2071, "qh_rename_adjacentvertex: merge #%d rename v%d (%d neighbors) to v%d (%d neighbors) dist %2.2g\n", zzval_(Ztotmerge), oldvertex->id, oldsize, newvertex->id, newsize, dist); } if (oldvertex->deleted || newvertex->deleted) { if (istrace || qh IStracing >= 2) { qh_fprintf(qh ferr, 2072, "qh_rename_adjacentvertex: ignore rename. Either v%d (%d) or v%d (%d) is deleted\n", oldvertex->id, oldvertex->deleted, newvertex->id, newvertex->deleted); } return; } if (oldsize == 0 || newsize == 0) { qh_fprintf(qh ferr, 2072, "qhull internal error (qh_rename_adjacentvertex): expecting neighbor facets for v%d and v%d. Got %d and %d neighbors resp.\n", oldvertex->id, newvertex->id, oldsize, newsize); qh_errexit(qh_ERRqhull, NULL, NULL); } FOREACHneighbor_(oldvertex) { if (neighbor->simplicial) { if (qh_setin(neighbor->vertices, newvertex)) { if (istrace || qh IStracing >= 2) { qh_fprintf(qh ferr, 2070, "qh_rename_adjacentvertex: simplicial f%d contains old v%d and new v%d. Will be marked degenerate by qh_renamevertex\n", neighbor->id, oldvertex->id, newvertex->id); } qh_makeridges(neighbor); /* no longer simplicial, nummerge==0, skipped by qh_maybe_duplicateridge */ }else { qh_replacefacetvertex(neighbor, oldvertex, newvertex); qh_setunique(&newvertex->neighbors, neighbor); qh_newvertices(neighbor->vertices); /* for qh_update_vertexneighbors of vertex neighbors */ } } } ridges= qh_vertexridges(oldvertex, qh_ALL); if (istrace) { FOREACHridge_(ridges) { qh_printridge(qh ferr, ridge); } } FOREACHneighbor_(oldvertex) { if (!neighbor->simplicial){ qh_addfacetvertex(neighbor, newvertex); qh_setunique(&newvertex->neighbors, neighbor); qh_newvertices(neighbor->vertices); /* for qh_update_vertexneighbors of vertex neighbors */ if (qh newfacet_list == qh facet_tail) { qh_removefacet(neighbor); /* add a neighbor to newfacet_list so that qh_partitionvisible has a newfacet */ qh_appendfacet(neighbor); neighbor->newfacet= True; } } } qh_renamevertex(oldvertex, newvertex, ridges, NULL, NULL); /* ridges invalidated */ if (oldvertex->deleted && !oldvertex->partitioned) { FOREACHneighbor_(newvertex) { if (!neighbor->visible) { qh_distplane(oldvertex->point, neighbor, &dist2); if (dist2>maxdist2) { maxdist2= dist2; maxfacet= neighbor; } } } trace2((qh ferr, 2096, "qh_rename_adjacentvertex: partition old p%d(v%d) as a coplanar point for furthest f%d dist %2.2g. Maybe repartition later (QH0031)\n", qh_pointid(oldvertex->point), oldvertex->id, maxfacet->id, maxdist2)) qh_partitioncoplanar(oldvertex->point, maxfacet, NULL, !qh_ALL); /* faster with maxdist2, otherwise duplicates distance tests from maxdist2/dist2 */ oldvertex->partitioned= True; } qh_settempfree(&ridges); } /* rename_adjacentvertex */ /*--------------------------------- qh_rename_sharedvertex( vertex, facet ) detect and rename if shared vertex in facet vertices have full ->neighbors returns: newvertex or NULL the vertex may still exist in other facets (i.e., a neighbor was pinched) does not change facet->neighbors updates vertex->neighbors notes: only called by qh_reducevertices after qh_remove_extravertices so f.vertices does not contain extraneous vertices a shared vertex for a facet is only in ridges to one neighbor this may undo a pinched facet it does not catch pinches involving multiple facets. These appear to be difficult to detect, since an exhaustive search is too expensive. design: if vertex only has two neighbors determine the ridges that contain the vertex determine the vertices shared by both neighbors if can find a new vertex in this set rename the vertex to the new vertex */ vertexT *qh_rename_sharedvertex(vertexT *vertex, facetT *facet) { facetT *neighbor, **neighborp, *neighborA= NULL; setT *vertices, *ridges; vertexT *newvertex= NULL; if (qh_setsize(vertex->neighbors) == 2) { neighborA= SETfirstt_(vertex->neighbors, facetT); if (neighborA == facet) neighborA= SETsecondt_(vertex->neighbors, facetT); }else if (qh hull_dim == 3) return NULL; else { qh visit_id++; FOREACHneighbor_(facet) neighbor->visitid= qh visit_id; FOREACHneighbor_(vertex) { if (neighbor->visitid == qh visit_id) { if (neighborA) return NULL; neighborA= neighbor; } } } if (!neighborA) { qh_fprintf(qh ferr, 6101, "qhull internal error (qh_rename_sharedvertex): v%d's neighbors not in f%d\n", vertex->id, facet->id); qh_errprint("ERRONEOUS", facet, NULL, NULL, vertex); qh_errexit(qh_ERRqhull, NULL, NULL); } if (neighborA) { /* avoid warning */ /* the vertex is shared by facet and neighborA */ ridges= qh_settemp(qh TEMPsize); neighborA->visitid= ++qh visit_id; qh_vertexridges_facet(vertex, facet, &ridges); trace2((qh ferr, 2037, "qh_rename_sharedvertex: p%d(v%d) is shared by f%d(%d ridges) and f%d\n", qh_pointid(vertex->point), vertex->id, facet->id, qh_setsize(ridges), neighborA->id)); zinc_(Zintersectnum); vertices= qh_vertexintersect_new(facet->vertices, neighborA->vertices); qh_setdel(vertices, vertex); qh_settemppush(vertices); if ((newvertex= qh_find_newvertex(vertex, vertices, ridges))) qh_renamevertex(vertex, newvertex, ridges, facet, neighborA); /* ridges invalidated */ qh_settempfree(&vertices); qh_settempfree(&ridges); } return newvertex; } /* rename_sharedvertex */ /*--------------------------------- qh_renameridgevertex( ridge, oldvertex, newvertex ) renames oldvertex as newvertex in ridge returns: True if renames oldvertex False if deleted the ridge notes: called by qh_renamevertex caller sets newvertex->delridge for deleted ridges (qh_reducevertices) design: delete oldvertex from ridge if newvertex already in ridge copy ridge->noconvex to another ridge if possible delete the ridge else insert newvertex into the ridge adjust the ridge's orientation */ boolT qh_renameridgevertex(ridgeT *ridge, vertexT *oldvertex, vertexT *newvertex) { int nth= 0, oldnth; facetT *temp; vertexT *vertex, **vertexp; oldnth= qh_setindex(ridge->vertices, oldvertex); if (oldnth < 0) { qh_fprintf(qh ferr, 6424, "qhull internal error (qh_renameridgevertex): oldvertex v%d not found in r%d. Cannot rename to v%d\n", oldvertex->id, ridge->id, newvertex->id); qh_errexit(qh_ERRqhull, NULL, ridge); } qh_setdelnthsorted(ridge->vertices, oldnth); FOREACHvertex_(ridge->vertices) { if (vertex == newvertex) { zinc_(Zdelridge); if (ridge->nonconvex) /* only one ridge has nonconvex set */ qh_copynonconvex(ridge); trace2((qh ferr, 2038, "qh_renameridgevertex: ridge r%d deleted. It contained both v%d and v%d\n", ridge->id, oldvertex->id, newvertex->id)); qh_delridge_merge(ridge); /* ridge.vertices deleted */ return False; } if (vertex->id < newvertex->id) break; nth++; } qh_setaddnth(&ridge->vertices, nth, newvertex); ridge->simplicialtop= False; ridge->simplicialbot= False; if (abs(oldnth - nth)%2) { trace3((qh ferr, 3010, "qh_renameridgevertex: swapped the top and bottom of ridge r%d\n", ridge->id)); temp= ridge->top; ridge->top= ridge->bottom; ridge->bottom= temp; } return True; } /* renameridgevertex */ /*--------------------------------- qh_renamevertex( oldvertex, newvertex, ridges, oldfacet, neighborA ) renames oldvertex as newvertex in ridges of non-simplicial neighbors set oldfacet/neighborA if oldvertex is shared between two facets (qh_rename_sharedvertex) otherwise qh_redundant_vertex or qh_rename_adjacentvertex returns: if oldfacet and multiple neighbors, oldvertex may still exist afterwards otherwise sets oldvertex->deleted for later deletion one or more ridges maybe deleted ridges is invalidated merges may be added to degen_mergeset via qh_maydropneighbor or qh_degen_redundant_facet notes: qh_rename_sharedvertex can not change neighbors of newvertex (since it's a subset) qh_redundant_vertex due to vertex->delridge for qh_reducevertices qh_rename_adjacentvertex for complete renames design: for each ridge in ridges rename oldvertex to newvertex and delete degenerate ridges if oldfacet not defined for each non-simplicial neighbor of oldvertex delete oldvertex from neighbor's vertices remove extra vertices from neighbor add oldvertex to qh.del_vertices else if oldvertex only between oldfacet and neighborA delete oldvertex from oldfacet and neighborA add oldvertex to qh.del_vertices else oldvertex is in oldfacet and neighborA and other facets (i.e., pinched) delete oldvertex from oldfacet delete oldfacet from old vertex's neighbors remove extra vertices (e.g., oldvertex) from neighborA */ void qh_renamevertex(vertexT *oldvertex, vertexT *newvertex, setT *ridges, facetT *oldfacet, facetT *neighborA) { facetT *neighbor, **neighborp; ridgeT *ridge, **ridgep; int topsize, bottomsize; boolT istrace= False; #ifndef qh_NOtrace if (qh IStracing >= 2 || oldvertex->id == qh tracevertex_id || newvertex->id == qh tracevertex_id) { istrace= True; qh_fprintf(qh ferr, 2086, "qh_renamevertex: rename v%d to v%d in %d ridges with old f%d and neighbor f%d\n", oldvertex->id, newvertex->id, qh_setsize(ridges), getid_(oldfacet), getid_(neighborA)); } #endif FOREACHridge_(ridges) { if (qh_renameridgevertex(ridge, oldvertex, newvertex)) { /* ridge is deleted if False, invalidating ridges */ topsize= qh_setsize(ridge->top->vertices); bottomsize= qh_setsize(ridge->bottom->vertices); if (topsize < qh hull_dim || (topsize == qh hull_dim && !ridge->top->simplicial && qh_setin(ridge->top->vertices, newvertex))) { trace4((qh ferr, 4070, "qh_renamevertex: ignore duplicate check for r%d. top f%d (size %d) will be degenerate after rename v%d to v%d\n", ridge->id, ridge->top->id, topsize, oldvertex->id, newvertex->id)); }else if (bottomsize < qh hull_dim || (bottomsize == qh hull_dim && !ridge->bottom->simplicial && qh_setin(ridge->bottom->vertices, newvertex))) { trace4((qh ferr, 4071, "qh_renamevertex: ignore duplicate check for r%d. bottom f%d (size %d) will be degenerate after rename v%d to v%d\n", ridge->id, ridge->bottom->id, bottomsize, oldvertex->id, newvertex->id)); }else qh_maybe_duplicateridge(ridge); } } if (!oldfacet) { /* stat Zrenameall or Zpinchduplicate */ if (istrace) qh_fprintf(qh ferr, 2087, "qh_renamevertex: renaming v%d to v%d in several facets for qh_redundant_vertex or MRGsubridge\n", oldvertex->id, newvertex->id); FOREACHneighbor_(oldvertex) { if (neighbor->simplicial) { qh_degen_redundant_facet(neighbor); /* e.g., rbox 175 C3,2e-13 D4 t1545235541 | qhull d */ }else { if (istrace) qh_fprintf(qh ferr, 4080, "qh_renamevertex: rename vertices in non-simplicial neighbor f%d of v%d\n", neighbor->id, oldvertex->id); qh_maydropneighbor(neighbor); qh_setdelsorted(neighbor->vertices, oldvertex); /* if degenerate, qh_degen_redundant_facet will add to mergeset */ if (qh_remove_extravertices(neighbor)) neighborp--; /* neighbor deleted from oldvertex neighborset */ qh_degen_redundant_facet(neighbor); /* either direction may be redundant, faster if combine? */ qh_test_redundant_neighbors(neighbor); qh_test_degen_neighbors(neighbor); } } if (!oldvertex->deleted) { oldvertex->deleted= True; qh_setappend(&qh del_vertices, oldvertex); } }else if (qh_setsize(oldvertex->neighbors) == 2) { zinc_(Zrenameshare); if (istrace) qh_fprintf(qh ferr, 3039, "qh_renamevertex: renaming v%d to v%d in oldfacet f%d for qh_rename_sharedvertex\n", oldvertex->id, newvertex->id, oldfacet->id); FOREACHneighbor_(oldvertex) { qh_setdelsorted(neighbor->vertices, oldvertex); qh_degen_redundant_facet(neighbor); } oldvertex->deleted= True; qh_setappend(&qh del_vertices, oldvertex); }else { zinc_(Zrenamepinch); if (istrace || qh IStracing >= 1) qh_fprintf(qh ferr, 3040, "qh_renamevertex: renaming pinched v%d to v%d between f%d and f%d\n", oldvertex->id, newvertex->id, oldfacet->id, neighborA->id); qh_setdelsorted(oldfacet->vertices, oldvertex); qh_setdel(oldvertex->neighbors, oldfacet); if (qh_remove_extravertices(neighborA)) qh_degen_redundant_facet(neighborA); } if (oldfacet) qh_degen_redundant_facet(oldfacet); } /* renamevertex */ /*--------------------------------- qh_test_appendmerge( facet, neighbor, simplicial ) test convexity and append to qh.facet_mergeset if non-convex if pre-merging, no-op if qh.SKIPconvex, or qh.MERGEexact and coplanar if simplicial, assumes centrum test is valid (e.g., adjacent, simplicial new facets) returns: true if appends facet/neighbor to qh.facet_mergeset sets facet->center as needed does not change facet->seen notes: called from qh_getmergeset_initial, qh_getmergeset, and qh_test_vneighbors must be at least as strong as qh_checkconvex (poly2.c) assumes !f.flipped design: exit if qh.SKIPconvex ('Q0') and !qh.POSTmerging if qh.cos_max ('An') is defined and merging coplanars if the angle between facet normals is too shallow append an angle-coplanar merge to qh.mergeset return True test convexity of facet and neighbor */ boolT qh_test_appendmerge(facetT *facet, facetT *neighbor, boolT simplicial) { realT angle= -REALmax; boolT okangle= False; if (qh SKIPconvex && !qh POSTmerging) return False; if (qh cos_max < REALmax/2 && (!qh MERGEexact || qh POSTmerging)) { angle= qh_getangle(facet->normal, neighbor->normal); okangle= True; zinc_(Zangletests); if (angle > qh cos_max) { zinc_(Zcoplanarangle); qh_appendmergeset(facet, neighbor, MRGanglecoplanar, 0.0, angle); trace2((qh ferr, 2039, "qh_test_appendmerge: coplanar angle %4.4g between f%d and f%d\n", angle, facet->id, neighbor->id)); return True; } } if (simplicial || qh hull_dim <= 3) return qh_test_centrum_merge(facet, neighbor, angle, okangle); else return qh_test_nonsimplicial_merge(facet, neighbor, angle, okangle); } /* test_appendmerge */ /*--------------------------------- qh_test_centrum_merge( facet, neighbor, angle, okangle ) test centrum convexity and append non-convex facets to qh.facet_mergeset 'angle' is angle between facets if okangle is true, otherwise use 0.0 returns: true if append facet/neighbor to qh.facet_mergeset sets facet->center as needed does not change facet->seen notes: called from test_appendmerge if adjacent simplicial facets or 2-d/3-d at least as strict as qh_checkconvex, including qh.DISTround ('En' and 'Rn') design: make facet's centrum if needed if facet's centrum is above the neighbor (qh.centrum_radius) set isconcave if facet's centrum is not below the neighbor (-qh.centrum_radius) set iscoplanar make neighbor's centrum if needed if neighbor's centrum is above the facet set isconcave else if neighbor's centrum is not below the facet set iscoplanar if isconcave or iscoplanar and merging coplanars get angle if needed (qh.ANGLEmerge 'An') append concave-coplanar, concave ,or coplanar merge to qh.mergeset */ boolT qh_test_centrum_merge(facetT *facet, facetT *neighbor, realT angle, boolT okangle) { coordT dist, dist2, mergedist; boolT isconcave= False, iscoplanar= False; if (!facet->center) facet->center= qh_getcentrum(facet); zzinc_(Zcentrumtests); qh_distplane(facet->center, neighbor, &dist); if (dist > qh centrum_radius) isconcave= True; else if (dist >= -qh centrum_radius) iscoplanar= True; if (!neighbor->center) neighbor->center= qh_getcentrum(neighbor); zzinc_(Zcentrumtests); qh_distplane(neighbor->center, facet, &dist2); if (dist2 > qh centrum_radius) isconcave= True; else if (!iscoplanar && dist2 >= -qh centrum_radius) iscoplanar= True; if (!isconcave && (!iscoplanar || (qh MERGEexact && !qh POSTmerging))) return False; if (!okangle && qh ANGLEmerge) { angle= qh_getangle(facet->normal, neighbor->normal); zinc_(Zangletests); } if (isconcave && iscoplanar) { zinc_(Zconcavecoplanarridge); if (dist > dist2) qh_appendmergeset(facet, neighbor, MRGconcavecoplanar, dist, angle); else qh_appendmergeset(neighbor, facet, MRGconcavecoplanar, dist2, angle); trace0((qh ferr, 36, "qh_test_centrum_merge: concave f%d to coplanar f%d, dist %4.4g and reverse dist %4.4g, angle %4.4g during p%d\n", facet->id, neighbor->id, dist, dist2, angle, qh furthest_id)); }else if (isconcave) { mergedist= fmax_(dist, dist2); zinc_(Zconcaveridge); qh_appendmergeset(facet, neighbor, MRGconcave, mergedist, angle); trace0((qh ferr, 37, "qh_test_centrum_merge: concave f%d to f%d, dist %4.4g and reverse dist %4.4g, angle %4.4g during p%d\n", facet->id, neighbor->id, dist, dist2, angle, qh furthest_id)); }else /* iscoplanar */ { mergedist= fmin_(fabs_(dist), fabs_(dist2)); zinc_(Zcoplanarcentrum); qh_appendmergeset(facet, neighbor, MRGcoplanar, mergedist, angle); trace2((qh ferr, 2097, "qh_test_centrum_merge: coplanar f%d to f%d dist %4.4g, reverse dist %4.4g angle %4.4g\n", facet->id, neighbor->id, dist, dist2, angle)); } return True; } /* test_centrum_merge */ /*--------------------------------- qh_test_degen_neighbors( facet ) append degenerate neighbors to qh.degen_mergeset notes: called at end of qh_mergefacet() and qh_renamevertex() call after test_redundant_facet() since MRGredundant is less expensive then MRGdegen a degenerate facet has fewer than hull_dim neighbors see: qh_merge_degenredundant() */ void qh_test_degen_neighbors(facetT *facet) { facetT *neighbor, **neighborp; int size; trace4((qh ferr, 4073, "qh_test_degen_neighbors: test for degenerate neighbors of f%d\n", facet->id)); FOREACHneighbor_(facet) { if (neighbor->visible) { qh_fprintf(qh ferr, 6359, "qhull internal error (qh_test_degen_neighbors): facet f%d has deleted neighbor f%d (qh.visible_list)\n", facet->id, neighbor->id); qh_errexit2(qh_ERRqhull, facet, neighbor); } if (neighbor->degenerate || neighbor->redundant || neighbor->dupridge) /* will merge or delete */ continue; /* merge flipped-degenerate facet before flipped facets */ if ((size= qh_setsize(neighbor->neighbors)) < qh hull_dim) { qh_appendmergeset(neighbor, neighbor, MRGdegen, 0.0, 1.0); trace2((qh ferr, 2019, "qh_test_degen_neighbors: f%d is degenerate with %d neighbors. Neighbor of f%d.\n", neighbor->id, size, facet->id)); } } } /* test_degen_neighbors */ /*--------------------------------- qh_test_nonsimplicial_merge( facet, neighbor, angle, okangle ) test centrum and vertex convexity and append non-convex or redundant facets to qh.facet_mergeset 'angle' is angle between facets if okangle is true, otherwise use 0.0 skips coplanar merges if pre-merging with qh.MERGEexact ('Qx') returns: true if appends facet/neighbor to qh.facet_mergeset sets facet->center as needed does not change facet->seen notes: only called from test_appendmerge if a non-simplicial facet and at least 4-d at least as strict as qh_checkconvex, including qh.DISTround ('En' and 'Rn') centrums must be < -qh.centrum_radius tests vertices as well as centrums since a facet may be twisted relative to its neighbor design: set precision constants for maxoutside, clearlyconcave, minvertex, and coplanarcentrum use maxoutside for coplanarcentrum if premerging with 'Qx' and qh_MAXcoplanarcentrum merges otherwise use qh.centrum_radious for coplanarcentrum make facet and neighbor centrums if needed isconcave if a centrum is above neighbor (coplanarcentrum) iscoplanar if a centrum is not below neighbor (-qh.centrum_radius) maybeconvex if a centrum is clearly below neighbor (-clearyconvex) return False if both centrums clearly below neighbor (-clearyconvex) return MRGconcave if isconcave facets are neither clearly convex nor clearly concave test vertices as well as centrums if maybeconvex determine mindist and maxdist for vertices of the other facet maybe MRGredundant otherwise determine mindist and maxdist for vertices of either facet maybe MRGredundant maybeconvex if a vertex is clearly below neighbor (-clearconvex) vertices are concave if dist > clearlyconcave vertices are twisted if dist > maxoutside (isconcave and maybeconvex) return False if not concave and pre-merge of 'Qx' (qh.MERGEexact) vertices are coplanar if dist in -minvertex..maxoutside if !isconcave, vertices are coplanar if dist >= -qh.MAXcoplanar (n*qh.premerge_centrum) return False if neither concave nor coplanar return MRGtwisted if isconcave and maybeconvex return MRGconcavecoplanar if isconcave and isconvex return MRGconcave if isconcave return MRGcoplanar if iscoplanar */ boolT qh_test_nonsimplicial_merge(facetT *facet, facetT *neighbor, realT angle, boolT okangle) { coordT dist, mindist, maxdist, mindist2, maxdist2, dist2, maxoutside, clearlyconcave, minvertex, clearlyconvex, mergedist, coplanarcentrum; boolT isconcave= False, iscoplanar= False, maybeconvex= False, isredundant= False; vertexT *maxvertex= NULL, *maxvertex2= NULL; maxoutside= fmax_(neighbor->maxoutside, qh ONEmerge + qh DISTround); maxoutside= fmax_(maxoutside, facet->maxoutside); clearlyconcave= qh_RATIOconcavehorizon * maxoutside; minvertex= fmax_(-qh min_vertex, qh MAXcoplanar); /* non-negative, not available per facet, not used for iscoplanar */ clearlyconvex= qh_RATIOconvexmerge * minvertex; /* must be convex for MRGtwisted */ if (qh MERGEexact && !qh POSTmerging && (facet->nummerge > qh_MAXcoplanarcentrum || neighbor->nummerge > qh_MAXcoplanarcentrum)) coplanarcentrum= maxoutside; else coplanarcentrum= qh centrum_radius; if (!facet->center) facet->center= qh_getcentrum(facet); zzinc_(Zcentrumtests); qh_distplane(facet->center, neighbor, &dist); if (dist > coplanarcentrum) isconcave= True; else if (dist >= -qh centrum_radius) iscoplanar= True; else if (dist < -clearlyconvex) maybeconvex= True; if (!neighbor->center) neighbor->center= qh_getcentrum(neighbor); zzinc_(Zcentrumtests); qh_distplane(neighbor->center, facet, &dist2); if (dist2 > coplanarcentrum) isconcave= True; else if (dist2 >= -qh centrum_radius) iscoplanar= True; else if (dist2 < -clearlyconvex) { if (maybeconvex) return False; /* both centrums clearly convex */ maybeconvex= True; } if (isconcave) { if (!okangle && qh ANGLEmerge) { angle= qh_getangle(facet->normal, neighbor->normal); zinc_(Zangletests); } mergedist= fmax_(dist, dist2); zinc_(Zconcaveridge); qh_appendmergeset(facet, neighbor, MRGconcave, mergedist, angle); trace0((qh ferr, 18, "qh_test_nonsimplicial_merge: concave centrum for f%d or f%d, dist %4.4g and reverse dist %4.4g, angle %4.4g during p%d\n", facet->id, neighbor->id, dist, dist2, angle, qh furthest_id)); return True; } /* neither clearly convex nor clearly concave, test vertices as well as centrums */ if (maybeconvex) { if (dist < -clearlyconvex) { maxdist= dist; /* facet centrum clearly convex, no need to test its vertex distance */ mindist= dist; maxvertex2= qh_furthestvertex(neighbor, facet, &maxdist2, &mindist2); if (!maxvertex2) { qh_appendmergeset(neighbor, facet, MRGredundant, maxdist2, qh_ANGLEnone); isredundant= True; } }else { /* dist2 < -clearlyconvex */ maxdist2= dist2; /* neighbor centrum clearly convex, no need to test its vertex distance */ mindist2= dist2; maxvertex= qh_furthestvertex(facet, neighbor, &maxdist, &mindist); if (!maxvertex) { qh_appendmergeset(facet, neighbor, MRGredundant, maxdist, qh_ANGLEnone); isredundant= True; } } }else { maxvertex= qh_furthestvertex(facet, neighbor, &maxdist, &mindist); if (maxvertex) { maxvertex2= qh_furthestvertex(neighbor, facet, &maxdist2, &mindist2); if (!maxvertex2) { qh_appendmergeset(neighbor, facet, MRGredundant, maxdist2, qh_ANGLEnone); isredundant= True; }else if (mindist < -clearlyconvex || mindist2 < -clearlyconvex) maybeconvex= True; }else { /* !maxvertex */ qh_appendmergeset(facet, neighbor, MRGredundant, maxdist, qh_ANGLEnone); isredundant= True; } } if (isredundant) { zinc_(Zredundantmerge); return True; } if (maxdist > clearlyconcave || maxdist2 > clearlyconcave) isconcave= True; else if (maybeconvex) { if (maxdist > maxoutside || maxdist2 > maxoutside) isconcave= True; /* MRGtwisted */ } if (!isconcave && qh MERGEexact && !qh POSTmerging) return False; if (isconcave && !iscoplanar) { if (maxdist < maxoutside && (-qh MAXcoplanar || (maxdist2 < maxoutside && mindist2 >= -qh MAXcoplanar))) iscoplanar= True; /* MRGconcavecoplanar */ }else if (!iscoplanar) { if (mindist >= -qh MAXcoplanar || mindist2 >= -qh MAXcoplanar) iscoplanar= True; /* MRGcoplanar */ } if (!isconcave && !iscoplanar) return False; if (!okangle && qh ANGLEmerge) { angle= qh_getangle(facet->normal, neighbor->normal); zinc_(Zangletests); } if (isconcave && maybeconvex) { zinc_(Ztwistedridge); if (maxdist > maxdist2) qh_appendmergeset(facet, neighbor, MRGtwisted, maxdist, angle); else qh_appendmergeset(neighbor, facet, MRGtwisted, maxdist2, angle); trace0((qh ferr, 27, "qh_test_nonsimplicial_merge: twisted concave f%d v%d to f%d v%d, dist %4.4g and reverse dist %4.4g, angle %4.4g during p%d\n", facet->id, getid_(maxvertex), neighbor->id, getid_(maxvertex2), maxdist, maxdist2, angle, qh furthest_id)); }else if (isconcave && iscoplanar) { zinc_(Zconcavecoplanarridge); if (maxdist > maxdist2) qh_appendmergeset(facet, neighbor, MRGconcavecoplanar, maxdist, angle); else qh_appendmergeset(neighbor, facet, MRGconcavecoplanar, maxdist2, angle); trace0((qh ferr, 28, "qh_test_nonsimplicial_merge: concave coplanar f%d v%d to f%d v%d, dist %4.4g and reverse dist %4.4g, angle %4.4g during p%d\n", facet->id, getid_(maxvertex), neighbor->id, getid_(maxvertex2), maxdist, maxdist2, angle, qh furthest_id)); }else if (isconcave) { mergedist= fmax_(maxdist, maxdist2); zinc_(Zconcaveridge); qh_appendmergeset(facet, neighbor, MRGconcave, mergedist, angle); trace0((qh ferr, 29, "qh_test_nonsimplicial_merge: concave f%d v%d to f%d v%d, dist %4.4g and reverse dist %4.4g, angle %4.4g during p%d\n", facet->id, getid_(maxvertex), neighbor->id, getid_(maxvertex2), maxdist, maxdist2, angle, qh furthest_id)); }else /* iscoplanar */ { mergedist= fmax_(fmax_(maxdist, maxdist2), fmax_(-mindist, -mindist2)); zinc_(Zcoplanarcentrum); qh_appendmergeset(facet, neighbor, MRGcoplanar, mergedist, angle); trace2((qh ferr, 2099, "qh_test_nonsimplicial_merge: coplanar f%d v%d to f%d v%d, dist %4.4g and reverse dist %4.4g, angle %4.4g during p%d\n", facet->id, getid_(maxvertex), neighbor->id, getid_(maxvertex2), maxdist, maxdist2, angle, qh furthest_id)); } return True; } /* test_nonsimplicial_merge */ /*--------------------------------- qh_test_redundant_neighbors( facet ) append degenerate facet or its redundant neighbors to qh.degen_mergeset returns: bumps vertex_visit notes: called at end of qh_mergefacet(), qh_mergecycle_all(), and qh_renamevertex call before qh_test_degen_neighbors (MRGdegen are more likely to cause problems) a redundant neighbor's vertices is a subset of the facet's vertices with pinched and flipped facets, a redundant neighbor may have a wildly different normal see qh_merge_degenredundant() and qh_-_facet() design: if facet is degenerate appends facet to degen_mergeset else appends redundant neighbors of facet to degen_mergeset */ void qh_test_redundant_neighbors(facetT *facet) { vertexT *vertex, **vertexp; facetT *neighbor, **neighborp; int size; trace4((qh ferr, 4022, "qh_test_redundant_neighbors: test neighbors of f%d vertex_visit %d\n", facet->id, qh vertex_visit+1)); if ((size= qh_setsize(facet->neighbors)) < qh hull_dim) { qh_appendmergeset(facet, facet, MRGdegen, 0.0, 1.0); trace2((qh ferr, 2017, "qh_test_redundant_neighbors: f%d is degenerate with %d neighbors.\n", facet->id, size)); }else { qh vertex_visit++; FOREACHvertex_(facet->vertices) vertex->visitid= qh vertex_visit; FOREACHneighbor_(facet) { if (neighbor->visible) { qh_fprintf(qh ferr, 6360, "qhull internal error (qh_test_redundant_neighbors): facet f%d has deleted neighbor f%d (qh.visible_list)\n", facet->id, neighbor->id); qh_errexit2(qh_ERRqhull, facet, neighbor); } if (neighbor->degenerate || neighbor->redundant || neighbor->dupridge) /* will merge or delete */ continue; if (facet->flipped && !neighbor->flipped) /* do not merge non-flipped into flipped */ continue; /* merge redundant-flipped facet first */ /* uses early out instead of checking vertex count */ FOREACHvertex_(neighbor->vertices) { if (vertex->visitid != qh vertex_visit) break; } if (!vertex) { qh_appendmergeset(neighbor, facet, MRGredundant, 0.0, 1.0); trace2((qh ferr, 2018, "qh_test_redundant_neighbors: f%d is contained in f%d. merge\n", neighbor->id, facet->id)); } } } } /* test_redundant_neighbors */ /*--------------------------------- qh_test_vneighbors( ) test vertex neighbors for convexity tests all facets on qh.newfacet_list returns: true if non-convex vneighbors appended to qh.facet_mergeset initializes vertex neighbors if needed notes: called by qh_all_merges from qh_postmerge if qh.TESTvneighbors ('Qv') assumes all facet neighbors have been tested this can be expensive this does not guarantee that a centrum is below all facets but it is unlikely uses qh.visit_id design: build vertex neighbors if necessary for all new facets for all vertices for each unvisited facet neighbor of the vertex test new facet and neighbor for convexity */ boolT qh_test_vneighbors(void /* qh.newfacet_list */) { facetT *newfacet, *neighbor, **neighborp; vertexT *vertex, **vertexp; int nummerges= 0; trace1((qh ferr, 1015, "qh_test_vneighbors: testing vertex neighbors for convexity\n")); if (!qh VERTEXneighbors) qh_vertexneighbors(); FORALLnew_facets newfacet->seen= False; FORALLnew_facets { newfacet->seen= True; newfacet->visitid= qh visit_id++; FOREACHneighbor_(newfacet) newfacet->visitid= qh visit_id; FOREACHvertex_(newfacet->vertices) { FOREACHneighbor_(vertex) { if (neighbor->seen || neighbor->visitid == qh visit_id) continue; if (qh_test_appendmerge(newfacet, neighbor, False)) /* ignores optimization for simplicial ridges */ nummerges++; } } } zadd_(Ztestvneighbor, nummerges); trace1((qh ferr, 1016, "qh_test_vneighbors: found %d non-convex, vertex neighbors\n", nummerges)); return (nummerges > 0); } /* test_vneighbors */ /*--------------------------------- qh_tracemerge( facet1, facet2 ) print trace message after merge */ void qh_tracemerge(facetT *facet1, facetT *facet2, mergeType mergetype) { boolT waserror= False; const char *mergename; #ifndef qh_NOtrace if(mergetype > 0 && mergetype < sizeof(mergetypes)/sizeof(char *)) mergename= mergetypes[mergetype]; else mergename= mergetypes[MRGnone]; if (qh IStracing >= 4) qh_errprint("MERGED", facet2, NULL, NULL, NULL); if (facet2 == qh tracefacet || (qh tracevertex && qh tracevertex->newfacet)) { qh_fprintf(qh ferr, 8085, "qh_tracemerge: trace facet and vertex after merge of f%d into f%d type %d (%s), furthest p%d\n", facet1->id, facet2->id, mergetype, mergename, qh furthest_id); if (facet2 != qh tracefacet) qh_errprint("TRACE", qh tracefacet, (qh tracevertex && qh tracevertex->neighbors) ? SETfirstt_(qh tracevertex->neighbors, facetT) : NULL, NULL, qh tracevertex); } if (qh tracevertex) { if (qh tracevertex->deleted) qh_fprintf(qh ferr, 8086, "qh_tracemerge: trace vertex deleted at furthest p%d\n", qh furthest_id); else qh_checkvertex(qh tracevertex, qh_ALL, &waserror); } if (qh tracefacet && qh tracefacet->normal && !qh tracefacet->visible) qh_checkfacet(qh tracefacet, True /* newmerge */, &waserror); #endif /* !qh_NOtrace */ if (qh CHECKfrequently || qh IStracing >= 4) { /* can't check polygon here */ if (qh IStracing >= 4 && qh num_facets < 500) { qh_printlists(); } qh_checkfacet(facet2, True /* newmerge */, &waserror); } if (waserror) qh_errexit(qh_ERRqhull, NULL, NULL); /* erroneous facet logged by qh_checkfacet */ } /* tracemerge */ /*--------------------------------- qh_tracemerging( ) print trace message during POSTmerging returns: updates qh.mergereport notes: called from qh_mergecycle() and qh_mergefacet() see: qh_buildtracing() */ void qh_tracemerging(void) { realT cpu; int total; time_t timedata; struct tm *tp; qh mergereport= zzval_(Ztotmerge); time(&timedata); tp= localtime(&timedata); cpu= qh_CPUclock; cpu /= qh_SECticks; total= zzval_(Ztotmerge) - zzval_(Zcyclehorizon) + zzval_(Zcyclefacettot); qh_fprintf(qh ferr, 8087, "\n\ At %d:%d:%d & %2.5g CPU secs, qhull has merged %d facets with max_outside %2.2g, min_vertex %2.2g.\n\ The hull contains %d facets and %d vertices.\n", tp->tm_hour, tp->tm_min, tp->tm_sec, cpu, total, qh max_outside, qh min_vertex, qh num_facets - qh num_visible, qh num_vertices-qh_setsize(qh del_vertices)); } /* tracemerging */ /*--------------------------------- qh_updatetested( facet1, facet2 ) clear facet2->tested and facet1->ridge->tested for merge returns: deletes facet2->center unless it's already large if so, clears facet2->ridge->tested notes: only called by qh_mergefacet design: clear facet2->tested clear ridge->tested for facet1's ridges if facet2 has a centrum if facet2 is large set facet2->keepcentrum else if facet2 has 3 vertices due to many merges, or not large and post merging clear facet2->keepcentrum unless facet2->keepcentrum clear facet2->center to recompute centrum later clear ridge->tested for facet2's ridges */ void qh_updatetested(facetT *facet1, facetT *facet2) { ridgeT *ridge, **ridgep; int size; facet2->tested= False; FOREACHridge_(facet1->ridges) ridge->tested= False; if (!facet2->center) return; size= qh_setsize(facet2->vertices); if (!facet2->keepcentrum) { if (size > qh hull_dim + qh_MAXnewcentrum) { facet2->keepcentrum= True; zinc_(Zwidevertices); } }else if (size <= qh hull_dim + qh_MAXnewcentrum) { /* center and keepcentrum was set */ if (size == qh hull_dim || qh POSTmerging) facet2->keepcentrum= False; /* if many merges need to recompute centrum */ } if (!facet2->keepcentrum) { qh_memfree(facet2->center, qh normal_size); facet2->center= NULL; FOREACHridge_(facet2->ridges) ridge->tested= False; } } /* updatetested */ /*--------------------------------- qh_vertexridges( vertex, allneighbors ) return temporary set of ridges adjacent to a vertex vertex->neighbors defined (qh_vertexneighbors) notes: uses qh.visit_id does not include implicit ridges for simplicial facets skips last neighbor, unless allneighbors. For new facets, the last neighbor shares ridges with adjacent neighbors if the last neighbor is not simplicial, it will have ridges for its simplicial neighbors Use allneighbors when a new cone is attached to an existing convex hull similar to qh_neighbor_vertices design: for each neighbor of vertex add ridges that include the vertex to ridges */ setT *qh_vertexridges(vertexT *vertex, boolT allneighbors) { facetT *neighbor, **neighborp; setT *ridges= qh_settemp(qh TEMPsize); int size; qh visit_id += 2; /* visit_id for vertex neighbors, visit_id-1 for facets of visited ridges */ FOREACHneighbor_(vertex) neighbor->visitid= qh visit_id; FOREACHneighbor_(vertex) { if (*neighborp || allneighbors) /* no new ridges in last neighbor */ qh_vertexridges_facet(vertex, neighbor, &ridges); } if (qh PRINTstatistics || qh IStracing) { size= qh_setsize(ridges); zinc_(Zvertexridge); zadd_(Zvertexridgetot, size); zmax_(Zvertexridgemax, size); trace3((qh ferr, 3011, "qh_vertexridges: found %d ridges for v%d\n", size, vertex->id)); } return ridges; } /* vertexridges */ /*--------------------------------- qh_vertexridges_facet( vertex, facet, ridges ) add adjacent ridges for vertex in facet neighbor->visitid==qh.visit_id if it hasn't been visited returns: ridges updated sets facet->visitid to qh.visit_id-1 design: for each ridge of facet if ridge of visited neighbor (i.e., unprocessed) if vertex in ridge append ridge mark facet processed */ void qh_vertexridges_facet(vertexT *vertex, facetT *facet, setT **ridges) { ridgeT *ridge, **ridgep; facetT *neighbor; int last_i= qh hull_dim-2; vertexT *second, *last; FOREACHridge_(facet->ridges) { neighbor= otherfacet_(ridge, facet); if (neighbor->visitid == qh visit_id) { if (SETfirst_(ridge->vertices) == vertex) { qh_setappend(ridges, ridge); }else if (last_i > 2) { second= SETsecondt_(ridge->vertices, vertexT); last= SETelemt_(ridge->vertices, last_i, vertexT); if (second->id >= vertex->id && last->id <= vertex->id) { /* vertices inverse sorted by id */ if (second == vertex || last == vertex) qh_setappend(ridges, ridge); else if (qh_setin(ridge->vertices, vertex)) qh_setappend(ridges, ridge); } }else if (SETelem_(ridge->vertices, last_i) == vertex || (last_i > 1 && SETsecond_(ridge->vertices) == vertex)) { qh_setappend(ridges, ridge); } } } facet->visitid= qh visit_id-1; } /* vertexridges_facet */ /*--------------------------------- qh_willdelete( facet, replace ) moves facet to visible list for qh_deletevisible sets facet->f.replace to replace (may be NULL) clears f.ridges and f.neighbors -- no longer valid returns: bumps qh.num_visible */ void qh_willdelete(facetT *facet, facetT *replace) { trace4((qh ferr, 4081, "qh_willdelete: move f%d to visible list, set its replacement as f%d, and clear f.neighbors and f.ridges\n", facet->id, getid_(replace))); if (!qh visible_list && qh newfacet_list) { qh_fprintf(qh ferr, 6378, "qhull internal error (qh_willdelete): expecting qh.visible_list at before qh.newfacet_list f%d. Got NULL\n", qh newfacet_list->id); qh_errexit2(qh_ERRqhull, NULL, NULL); } qh_removefacet(facet); qh_prependfacet(facet, &qh visible_list); qh num_visible++; facet->visible= True; facet->f.replace= replace; if (facet->ridges) SETfirst_(facet->ridges)= NULL; if (facet->neighbors) SETfirst_(facet->neighbors)= NULL; } /* willdelete */ #else /* qh_NOmerge */ void qh_all_vertexmerges(int apexpointid, facetT *facet, facetT **retryfacet) { QHULL_UNUSED(apexpointid) QHULL_UNUSED(facet) QHULL_UNUSED(retryfacet) } void qh_premerge(int apexpointid, realT maxcentrum, realT maxangle) { QHULL_UNUSED(apexpointid) QHULL_UNUSED(maxcentrum) QHULL_UNUSED(maxangle) } void qh_postmerge(const char *reason, realT maxcentrum, realT maxangle, boolT vneighbors) { QHULL_UNUSED(reason) QHULL_UNUSED(maxcentrum) QHULL_UNUSED(maxangle) QHULL_UNUSED(vneighbors) } void qh_checkdelfacet(facetT *facet, setT *mergeset) { QHULL_UNUSED(facet) QHULL_UNUSED(mergeset) } void qh_checkdelridge(void /* qh.visible_facets, vertex_mergeset */) { } boolT qh_checkzero(boolT testall) { QHULL_UNUSED(testall) return True; } void qh_freemergesets(void) { } void qh_initmergesets(void) { } void qh_merge_pinchedvertices(int apexpointid /* qh.newfacet_list */) { QHULL_UNUSED(apexpointid) } #endif /* qh_NOmerge */