diff options
Diffstat (limited to 'extern/qhull/src/io.c')
| -rw-r--r-- | extern/qhull/src/io.c | 4089 |
1 files changed, 4089 insertions, 0 deletions
diff --git a/extern/qhull/src/io.c b/extern/qhull/src/io.c new file mode 100644 index 00000000000..9b0ccdd0b24 --- /dev/null +++ b/extern/qhull/src/io.c @@ -0,0 +1,4089 @@ +/*<html><pre> -<a href="qh-io.htm" + >-------------------------------</a><a name="TOP">-</a> + + io.c + Input/Output routines of qhull application + + see qh-io.htm and io.h + + see user.c for qh_errprint and qh_printfacetlist + + unix.c calls qh_readpoints and qh_produce_output + + unix.c and user.c are the only callers of io.c functions + This allows the user to avoid loading io.o from qhull.a + + copyright (c) 1993-2002 The Geometry Center +*/ + +#include "qhull_a.h" + +/*========= -prototypes for internal functions ========= */ + +static int qh_compare_facetarea(const void *p1, const void *p2); +static int qh_compare_facetmerge(const void *p1, const void *p2); +static int qh_compare_facetvisit(const void *p1, const void *p2); +int qh_compare_vertexpoint(const void *p1, const void *p2); /* not used */ + +/*========= -functions in alphabetical order after qh_produce_output() =====*/ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="produce_output">-</a> + + qh_produce_output() + prints out the result of qhull in desired format + if qh.GETarea + computes and prints area and volume + qh.PRINTout[] is an array of output formats + + notes: + prints output in qh.PRINTout order +*/ +void qh_produce_output(void) { + int i, tempsize= qh_setsize ((setT*)qhmem.tempstack), d_1; + + if (qh VORONOI) { + qh_clearcenters (qh_ASvoronoi); + qh_vertexneighbors(); + } + if (qh TRIangulate) { + qh_triangulate(); + if (qh VERIFYoutput && !qh CHECKfrequently) + qh_checkpolygon (qh facet_list); + } + qh_findgood_all (qh facet_list); + if (qh GETarea) + qh_getarea(qh facet_list); + if (qh KEEParea || qh KEEPmerge || qh KEEPminArea < REALmax/2) + qh_markkeep (qh facet_list); + if (qh PRINTsummary) + qh_printsummary(qh ferr); + else if (qh PRINTout[0] == qh_PRINTnone) + qh_printsummary(qh fout); + for (i= 0; i < qh_PRINTEND; i++) + qh_printfacets (qh fout, qh PRINTout[i], qh facet_list, NULL, !qh_ALL); + qh_allstatistics(); + if (qh PRINTprecision && !qh MERGING && (qh JOGGLEmax > REALmax/2 || qh RERUN)) + qh_printstats (qh ferr, qhstat precision, NULL); + if (qh VERIFYoutput && (zzval_(Zridge) > 0 || zzval_(Zridgemid) > 0)) + qh_printstats (qh ferr, qhstat vridges, NULL); + if (qh PRINTstatistics) { + qh_collectstatistics(); + qh_printstatistics(qh ferr, ""); + qh_memstatistics (qh ferr); + d_1= sizeof(setT) + (qh hull_dim - 1) * SETelemsize; + fprintf(qh ferr, "\ + size in bytes: merge %ld ridge %ld vertex %ld facet %ld\n\ + normal %d ridge vertices %d facet vertices or neighbors %ld\n", + sizeof(mergeT), sizeof(ridgeT), + sizeof(vertexT), sizeof(facetT), + qh normal_size, d_1, d_1 + SETelemsize); + } + if (qh_setsize ((setT*)qhmem.tempstack) != tempsize) { + fprintf (qh ferr, "qhull internal error (qh_produce_output): temporary sets not empty (%d)\n", + qh_setsize ((setT*)qhmem.tempstack)); + qh_errexit (qh_ERRqhull, NULL, NULL); + } +} /* produce_output */ + + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="dfacet">-</a> + + dfacet( id ) + print facet by id, for debugging + +*/ +void dfacet (unsigned id) { + facetT *facet; + + FORALLfacets { + if (facet->id == id) { + qh_printfacet (qh fout, facet); + break; + } + } +} /* dfacet */ + + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="dvertex">-</a> + + dvertex( id ) + print vertex by id, for debugging +*/ +void dvertex (unsigned id) { + vertexT *vertex; + + FORALLvertices { + if (vertex->id == id) { + qh_printvertex (qh fout, vertex); + break; + } + } +} /* dvertex */ + + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="compare_vertexpoint">-</a> + + qh_compare_vertexpoint( p1, p2 ) + used by qsort() to order vertices by point id +*/ +int qh_compare_vertexpoint(const void *p1, const void *p2) { + vertexT *a= *((vertexT **)p1), *b= *((vertexT **)p2); + + return ((qh_pointid(a->point) > qh_pointid(b->point)?1:-1)); +} /* compare_vertexpoint */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="compare_facetarea">-</a> + + qh_compare_facetarea( p1, p2 ) + used by qsort() to order facets by area +*/ +static int qh_compare_facetarea(const void *p1, const void *p2) { + facetT *a= *((facetT **)p1), *b= *((facetT **)p2); + + if (!a->isarea) + return -1; + if (!b->isarea) + return 1; + if (a->f.area > b->f.area) + return 1; + else if (a->f.area == b->f.area) + return 0; + return -1; +} /* compare_facetarea */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="compare_facetmerge">-</a> + + qh_compare_facetmerge( p1, p2 ) + used by qsort() to order facets by number of merges +*/ +static int qh_compare_facetmerge(const void *p1, const void *p2) { + facetT *a= *((facetT **)p1), *b= *((facetT **)p2); + + return (a->nummerge - b->nummerge); +} /* compare_facetvisit */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="compare_facetvisit">-</a> + + qh_compare_facetvisit( p1, p2 ) + used by qsort() to order facets by visit id or id +*/ +static int qh_compare_facetvisit(const void *p1, const void *p2) { + facetT *a= *((facetT **)p1), *b= *((facetT **)p2); + int i,j; + + if (!(i= a->visitid)) + i= - a->id; /* do not convert to int */ + if (!(j= b->visitid)) + j= - b->id; + return (i - j); +} /* compare_facetvisit */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="countfacets">-</a> + + qh_countfacets( facetlist, facets, printall, + numfacets, numsimplicial, totneighbors, numridges, numcoplanar, numtricoplanars ) + count good facets for printing and set visitid + if allfacets, ignores qh_skipfacet() + + notes: + qh_printsummary and qh_countfacets must match counts + + returns: + numfacets, numsimplicial, total neighbors, numridges, coplanars + each facet with ->visitid indicating 1-relative position + ->visitid==0 indicates not good + + notes + numfacets >= numsimplicial + if qh.NEWfacets, + does not count visible facets (matches qh_printafacet) + + design: + for all facets on facetlist and in facets set + unless facet is skipped or visible (i.e., will be deleted) + mark facet->visitid + update counts +*/ +void qh_countfacets (facetT *facetlist, setT *facets, boolT printall, + int *numfacetsp, int *numsimplicialp, int *totneighborsp, int *numridgesp, int *numcoplanarsp, int *numtricoplanarsp) { + facetT *facet, **facetp; + int numfacets= 0, numsimplicial= 0, numridges= 0, totneighbors= 0, numcoplanars= 0, numtricoplanars= 0; + + FORALLfacet_(facetlist) { + if ((facet->visible && qh NEWfacets) + || (!printall && qh_skipfacet(facet))) + facet->visitid= 0; + else { + facet->visitid= ++numfacets; + totneighbors += qh_setsize (facet->neighbors); + if (facet->simplicial) { + numsimplicial++; + if (facet->keepcentrum && facet->tricoplanar) + numtricoplanars++; + }else + numridges += qh_setsize (facet->ridges); + if (facet->coplanarset) + numcoplanars += qh_setsize (facet->coplanarset); + } + } + FOREACHfacet_(facets) { + if ((facet->visible && qh NEWfacets) + || (!printall && qh_skipfacet(facet))) + facet->visitid= 0; + else { + facet->visitid= ++numfacets; + totneighbors += qh_setsize (facet->neighbors); + if (facet->simplicial){ + numsimplicial++; + if (facet->keepcentrum && facet->tricoplanar) + numtricoplanars++; + }else + numridges += qh_setsize (facet->ridges); + if (facet->coplanarset) + numcoplanars += qh_setsize (facet->coplanarset); + } + } + qh visit_id += numfacets+1; + *numfacetsp= numfacets; + *numsimplicialp= numsimplicial; + *totneighborsp= totneighbors; + *numridgesp= numridges; + *numcoplanarsp= numcoplanars; + *numtricoplanarsp= numtricoplanars; +} /* countfacets */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="detvnorm">-</a> + + qh_detvnorm( vertex, vertexA, centers, offset ) + compute separating plane of the Voronoi diagram for a pair of input sites + centers= set of facets (i.e., Voronoi vertices) + facet->visitid= 0 iff vertex-at-infinity (i.e., unbounded) + + assumes: + qh_ASvoronoi and qh_vertexneighbors() already set + + returns: + norm + a pointer into qh.gm_matrix to qh.hull_dim-1 reals + copy the data before reusing qh.gm_matrix + offset + if 'QVn' + sign adjusted so that qh.GOODvertexp is inside + else + sign adjusted so that vertex is inside + + qh.gm_matrix= simplex of points from centers relative to first center + + notes: + in io.c so that code for 'v Tv' can be removed by removing io.c + returns pointer into qh.gm_matrix to avoid tracking of temporary memory + + design: + determine midpoint of input sites + build points as the set of Voronoi vertices + select a simplex from points (if necessary) + include midpoint if the Voronoi region is unbounded + relocate the first vertex of the simplex to the origin + compute the normalized hyperplane through the simplex + orient the hyperplane toward 'QVn' or 'vertex' + if 'Tv' or 'Ts' + if bounded + test that hyperplane is the perpendicular bisector of the input sites + test that Voronoi vertices not in the simplex are still on the hyperplane + free up temporary memory +*/ +pointT *qh_detvnorm (vertexT *vertex, vertexT *vertexA, setT *centers, realT *offsetp) { + facetT *facet, **facetp; + int i, k, pointid, pointidA, point_i, point_n; + setT *simplex= NULL; + pointT *point, **pointp, *point0, *midpoint, *normal, *inpoint; + coordT *coord, *gmcoord, *normalp; + setT *points= qh_settemp (qh TEMPsize); + boolT nearzero= False; + boolT unbounded= False; + int numcenters= 0; + int dim= qh hull_dim - 1; + realT dist, offset, angle, zero= 0.0; + + midpoint= qh gm_matrix + qh hull_dim * qh hull_dim; /* last row */ + for (k= 0; k < dim; k++) + midpoint[k]= (vertex->point[k] + vertexA->point[k])/2; + FOREACHfacet_(centers) { + numcenters++; + if (!facet->visitid) + unbounded= True; + else { + if (!facet->center) + facet->center= qh_facetcenter (facet->vertices); + qh_setappend (&points, facet->center); + } + } + if (numcenters > dim) { + simplex= qh_settemp (qh TEMPsize); + qh_setappend (&simplex, vertex->point); + if (unbounded) + qh_setappend (&simplex, midpoint); + qh_maxsimplex (dim, points, NULL, 0, &simplex); + qh_setdelnth (simplex, 0); + }else if (numcenters == dim) { + if (unbounded) + qh_setappend (&points, midpoint); + simplex= points; + }else { + fprintf(qh ferr, "qh_detvnorm: too few points (%d) to compute separating plane\n", numcenters); + qh_errexit (qh_ERRqhull, NULL, NULL); + } + i= 0; + gmcoord= qh gm_matrix; + point0= SETfirstt_(simplex, pointT); + FOREACHpoint_(simplex) { + if (qh IStracing >= 4) + qh_printmatrix(qh ferr, "qh_detvnorm: Voronoi vertex or midpoint", + &point, 1, dim); + if (point != point0) { + qh gm_row[i++]= gmcoord; + coord= point0; + for (k= dim; k--; ) + *(gmcoord++)= *point++ - *coord++; + } + } + qh gm_row[i]= gmcoord; /* does not overlap midpoint, may be used later for qh_areasimplex */ + normal= gmcoord; + qh_sethyperplane_gauss (dim, qh gm_row, point0, True, + normal, &offset, &nearzero); + if (qh GOODvertexp == vertexA->point) + inpoint= vertexA->point; + else + inpoint= vertex->point; + zinc_(Zdistio); + dist= qh_distnorm (dim, inpoint, normal, &offset); + if (dist > 0) { + offset= -offset; + normalp= normal; + for (k= dim; k--; ) { + *normalp= -(*normalp); + normalp++; + } + } + if (qh VERIFYoutput || qh PRINTstatistics) { + pointid= qh_pointid (vertex->point); + pointidA= qh_pointid (vertexA->point); + if (!unbounded) { + zinc_(Zdiststat); + dist= qh_distnorm (dim, midpoint, normal, &offset); + if (dist < 0) + dist= -dist; + zzinc_(Zridgemid); + wwmax_(Wridgemidmax, dist); + wwadd_(Wridgemid, dist); + trace4((qh ferr, "qh_detvnorm: points %d %d midpoint dist %2.2g\n", + pointid, pointidA, dist)); + for (k= 0; k < dim; k++) + midpoint[k]= vertexA->point[k] - vertex->point[k]; /* overwrites midpoint! */ + qh_normalize (midpoint, dim, False); + angle= qh_distnorm (dim, midpoint, normal, &zero); /* qh_detangle uses dim+1 */ + if (angle < 0.0) + angle= angle + 1.0; + else + angle= angle - 1.0; + if (angle < 0.0) + angle -= angle; + trace4((qh ferr, "qh_detvnorm: points %d %d angle %2.2g nearzero %d\n", + pointid, pointidA, angle, nearzero)); + if (nearzero) { + zzinc_(Zridge0); + wwmax_(Wridge0max, angle); + wwadd_(Wridge0, angle); + }else { + zzinc_(Zridgeok) + wwmax_(Wridgeokmax, angle); + wwadd_(Wridgeok, angle); + } + } + if (simplex != points) { + FOREACHpoint_i_(points) { + if (!qh_setin (simplex, point)) { + facet= SETelemt_(centers, point_i, facetT); + zinc_(Zdiststat); + dist= qh_distnorm (dim, point, normal, &offset); + if (dist < 0) + dist= -dist; + zzinc_(Zridge); + wwmax_(Wridgemax, dist); + wwadd_(Wridge, dist); + trace4((qh ferr, "qh_detvnorm: points %d %d Voronoi vertex %d dist %2.2g\n", + pointid, pointidA, facet->visitid, dist)); + } + } + } + } + *offsetp= offset; + if (simplex != points) + qh_settempfree (&simplex); + qh_settempfree (&points); + return normal; +} /* detvnorm */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="detvridge">-</a> + + qh_detvridge( vertexA ) + determine Voronoi ridge from 'seen' neighbors of vertexA + include one vertex-at-infinite if an !neighbor->visitid + + returns: + temporary set of centers (facets, i.e., Voronoi vertices) + sorted by center id +*/ +setT *qh_detvridge (vertexT *vertex) { + setT *centers= qh_settemp (qh TEMPsize); + setT *tricenters= qh_settemp (qh TEMPsize); + facetT *neighbor, **neighborp; + boolT firstinf= True; + + FOREACHneighbor_(vertex) { + if (neighbor->seen) { + if (neighbor->visitid) { + if (!neighbor->tricoplanar || qh_setunique (&tricenters, neighbor->center)) + qh_setappend (¢ers, neighbor); + }else if (firstinf) { + firstinf= False; + qh_setappend (¢ers, neighbor); + } + } + } + qsort (SETaddr_(centers, facetT), qh_setsize (centers), + sizeof (facetT *), qh_compare_facetvisit); + qh_settempfree (&tricenters); + return centers; +} /* detvridge */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="detvridge3">-</a> + + qh_detvridge3( atvertex, vertex ) + determine 3-d Voronoi ridge from 'seen' neighbors of atvertex and vertex + include one vertex-at-infinite for !neighbor->visitid + assumes all facet->seen2= True + + returns: + temporary set of centers (facets, i.e., Voronoi vertices) + listed in adjacency order (not oriented) + all facet->seen2= True + + design: + mark all neighbors of atvertex + for each adjacent neighbor of both atvertex and vertex + if neighbor selected + add neighbor to set of Voronoi vertices +*/ +setT *qh_detvridge3 (vertexT *atvertex, vertexT *vertex) { + setT *centers= qh_settemp (qh TEMPsize); + setT *tricenters= qh_settemp (qh TEMPsize); + facetT *neighbor, **neighborp, *facet= NULL; + boolT firstinf= True; + + FOREACHneighbor_(atvertex) + neighbor->seen2= False; + FOREACHneighbor_(vertex) { + if (!neighbor->seen2) { + facet= neighbor; + break; + } + } + while (facet) { + facet->seen2= True; + if (neighbor->seen) { + if (facet->visitid) { + if (!facet->tricoplanar || qh_setunique (&tricenters, facet->center)) + qh_setappend (¢ers, facet); + }else if (firstinf) { + firstinf= False; + qh_setappend (¢ers, facet); + } + } + FOREACHneighbor_(facet) { + if (!neighbor->seen2) { + if (qh_setin (vertex->neighbors, neighbor)) + break; + else + neighbor->seen2= True; + } + } + facet= neighbor; + } + if (qh CHECKfrequently) { + FOREACHneighbor_(vertex) { + if (!neighbor->seen2) { + fprintf (stderr, "qh_detvridge3: neigbors of vertex p%d are not connected at facet %d\n", + qh_pointid (vertex->point), neighbor->id); + qh_errexit (qh_ERRqhull, neighbor, NULL); + } + } + } + FOREACHneighbor_(atvertex) + neighbor->seen2= True; + qh_settempfree (&tricenters); + return centers; +} /* detvridge3 */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="eachvoronoi">-</a> + + qh_eachvoronoi( fp, printvridge, vertex, visitall, innerouter, inorder ) + if visitall, + visit all Voronoi ridges for vertex (i.e., an input site) + else + visit all unvisited Voronoi ridges for vertex + all vertex->seen= False if unvisited + assumes + all facet->seen= False + all facet->seen2= True (for qh_detvridge3) + all facet->visitid == 0 if vertex_at_infinity + == index of Voronoi vertex + >= qh.num_facets if ignored + innerouter: + qh_RIDGEall-- both inner (bounded) and outer (unbounded) ridges + qh_RIDGEinner- only inner + qh_RIDGEouter- only outer + + if inorder + orders vertices for 3-d Voronoi diagrams + + returns: + number of visited ridges (does not include previously visited ridges) + + if printvridge, + calls printvridge( fp, vertex, vertexA, centers) + fp== any pointer (assumes FILE*) + vertex,vertexA= pair of input sites that define a Voronoi ridge + centers= set of facets (i.e., Voronoi vertices) + ->visitid == index or 0 if vertex_at_infinity + ordered for 3-d Voronoi diagram + notes: + uses qh.vertex_visit + + see: + qh_eachvoronoi_all() + + design: + mark selected neighbors of atvertex + for each selected neighbor (either Voronoi vertex or vertex-at-infinity) + for each unvisited vertex + if atvertex and vertex share more than d-1 neighbors + bump totalcount + if printvridge defined + build the set of shared neighbors (i.e., Voronoi vertices) + call printvridge +*/ +int qh_eachvoronoi (FILE *fp, printvridgeT printvridge, vertexT *atvertex, boolT visitall, qh_RIDGE innerouter, boolT inorder) { + boolT unbounded; + int count; + facetT *neighbor, **neighborp, *neighborA, **neighborAp; + setT *centers; + setT *tricenters= qh_settemp (qh TEMPsize); + + vertexT *vertex, **vertexp; + boolT firstinf; + unsigned int numfacets= (unsigned int)qh num_facets; + int totridges= 0; + + qh vertex_visit++; + atvertex->seen= True; + if (visitall) { + FORALLvertices + vertex->seen= False; + } + FOREACHneighbor_(atvertex) { + if (neighbor->visitid < numfacets) + neighbor->seen= True; + } + FOREACHneighbor_(atvertex) { + if (neighbor->seen) { + FOREACHvertex_(neighbor->vertices) { + if (vertex->visitid != qh vertex_visit && !vertex->seen) { + vertex->visitid= qh vertex_visit; + count= 0; + firstinf= True; + qh_settruncate (tricenters, 0); + FOREACHneighborA_(vertex) { + if (neighborA->seen) { + if (neighborA->visitid) { + if (!neighborA->tricoplanar || qh_setunique (&tricenters, neighborA->center)) + count++; + }else if (firstinf) { + count++; + firstinf= False; + } + } + } + if (count >= qh hull_dim - 1) { /* e.g., 3 for 3-d Voronoi */ + if (firstinf) { + if (innerouter == qh_RIDGEouter) + continue; + unbounded= False; + }else { + if (innerouter == qh_RIDGEinner) + continue; + unbounded= True; + } + totridges++; + trace4((qh ferr, "qh_eachvoronoi: Voronoi ridge of %d vertices between sites %d and %d\n", + count, qh_pointid (atvertex->point), qh_pointid (vertex->point))); + if (printvridge) { + if (inorder && qh hull_dim == 3+1) /* 3-d Voronoi diagram */ + centers= qh_detvridge3 (atvertex, vertex); + else + centers= qh_detvridge (vertex); + (*printvridge) (fp, atvertex, vertex, centers, unbounded); + qh_settempfree (¢ers); + } + } + } + } + } + } + FOREACHneighbor_(atvertex) + neighbor->seen= False; + qh_settempfree (&tricenters); + return totridges; +} /* eachvoronoi */ + + +/*-<a href="qh-poly.htm#TOC" + >-------------------------------</a><a name="eachvoronoi_all">-</a> + + qh_eachvoronoi_all( fp, printvridge, isupper, innerouter, inorder ) + visit all Voronoi ridges + + innerouter: + see qh_eachvoronoi() + + if inorder + orders vertices for 3-d Voronoi diagrams + + returns + total number of ridges + + if isupper == facet->upperdelaunay (i.e., a Vornoi vertex) + facet->visitid= Voronoi vertex index (same as 'o' format) + else + facet->visitid= 0 + + if printvridge, + calls printvridge( fp, vertex, vertexA, centers) + [see qh_eachvoronoi] + + notes: + Not used for qhull.exe + same effect as qh_printvdiagram but ridges not sorted by point id +*/ +int qh_eachvoronoi_all (FILE *fp, printvridgeT printvridge, boolT isupper, qh_RIDGE innerouter, boolT inorder) { + facetT *facet; + vertexT *vertex; + int numcenters= 1; /* vertex 0 is vertex-at-infinity */ + int totridges= 0; + + qh_clearcenters (qh_ASvoronoi); + qh_vertexneighbors(); + maximize_(qh visit_id, (unsigned) qh num_facets); + FORALLfacets { + facet->visitid= 0; + facet->seen= False; + facet->seen2= True; + } + FORALLfacets { + if (facet->upperdelaunay == isupper) + facet->visitid= numcenters++; + } + FORALLvertices + vertex->seen= False; + FORALLvertices { + if (qh GOODvertex > 0 && qh_pointid(vertex->point)+1 != qh GOODvertex) + continue; + totridges += qh_eachvoronoi (fp, printvridge, vertex, + !qh_ALL, innerouter, inorder); + } + return totridges; +} /* eachvoronoi_all */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="facet2point">-</a> + + qh_facet2point( facet, point0, point1, mindist ) + return two projected temporary vertices for a 2-d facet + may be non-simplicial + + returns: + point0 and point1 oriented and projected to the facet + returns mindist (maximum distance below plane) +*/ +void qh_facet2point(facetT *facet, pointT **point0, pointT **point1, realT *mindist) { + vertexT *vertex0, *vertex1; + realT dist; + + if (facet->toporient ^ qh_ORIENTclock) { + vertex0= SETfirstt_(facet->vertices, vertexT); + vertex1= SETsecondt_(facet->vertices, vertexT); + }else { + vertex1= SETfirstt_(facet->vertices, vertexT); + vertex0= SETsecondt_(facet->vertices, vertexT); + } + zadd_(Zdistio, 2); + qh_distplane(vertex0->point, facet, &dist); + *mindist= dist; + *point0= qh_projectpoint(vertex0->point, facet, dist); + qh_distplane(vertex1->point, facet, &dist); + minimize_(*mindist, dist); + *point1= qh_projectpoint(vertex1->point, facet, dist); +} /* facet2point */ + + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="facetvertices">-</a> + + qh_facetvertices( facetlist, facets, allfacets ) + returns temporary set of vertices in a set and/or list of facets + if allfacets, ignores qh_skipfacet() + + returns: + vertices with qh.vertex_visit + + notes: + optimized for allfacets of facet_list + + design: + if allfacets of facet_list + create vertex set from vertex_list + else + for each selected facet in facets or facetlist + append unvisited vertices to vertex set +*/ +setT *qh_facetvertices (facetT *facetlist, setT *facets, boolT allfacets) { + setT *vertices; + facetT *facet, **facetp; + vertexT *vertex, **vertexp; + + qh vertex_visit++; + if (facetlist == qh facet_list && allfacets && !facets) { + vertices= qh_settemp (qh num_vertices); + FORALLvertices { + vertex->visitid= qh vertex_visit; + qh_setappend (&vertices, vertex); + } + }else { + vertices= qh_settemp (qh TEMPsize); + FORALLfacet_(facetlist) { + if (!allfacets && qh_skipfacet (facet)) + continue; + FOREACHvertex_(facet->vertices) { + if (vertex->visitid != qh vertex_visit) { + vertex->visitid= qh vertex_visit; + qh_setappend (&vertices, vertex); + } + } + } + } + FOREACHfacet_(facets) { + if (!allfacets && qh_skipfacet (facet)) + continue; + FOREACHvertex_(facet->vertices) { + if (vertex->visitid != qh vertex_visit) { + vertex->visitid= qh vertex_visit; + qh_setappend (&vertices, vertex); + } + } + } + return vertices; +} /* facetvertices */ + +/*-<a href="qh-geom.htm#TOC" + >-------------------------------</a><a name="geomplanes">-</a> + + qh_geomplanes( facet, outerplane, innerplane ) + return outer and inner planes for Geomview + qh.PRINTradius is size of vertices and points (includes qh.JOGGLEmax) + + notes: + assume precise calculations in io.c with roundoff covered by qh_GEOMepsilon +*/ +void qh_geomplanes (facetT *facet, realT *outerplane, realT *innerplane) { + realT radius; + + if (qh MERGING || qh JOGGLEmax < REALmax/2) { + qh_outerinner (facet, outerplane, innerplane); + radius= qh PRINTradius; + if (qh JOGGLEmax < REALmax/2) + radius -= qh JOGGLEmax * sqrt (qh hull_dim); /* already accounted for in qh_outerinner() */ + *outerplane += radius; + *innerplane -= radius; + if (qh PRINTcoplanar || qh PRINTspheres) { + *outerplane += qh MAXabs_coord * qh_GEOMepsilon; + *innerplane -= qh MAXabs_coord * qh_GEOMepsilon; + } + }else + *innerplane= *outerplane= 0; +} /* geomplanes */ + + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="markkeep">-</a> + + qh_markkeep( facetlist ) + mark good facets that meet qh.KEEParea, qh.KEEPmerge, and qh.KEEPminArea + ignores visible facets (not part of convex hull) + + returns: + may clear facet->good + recomputes qh.num_good + + design: + get set of good facets + if qh.KEEParea + sort facets by area + clear facet->good for all but n largest facets + if qh.KEEPmerge + sort facets by merge count + clear facet->good for all but n most merged facets + if qh.KEEPminarea + clear facet->good if area too small + update qh.num_good +*/ +void qh_markkeep (facetT *facetlist) { + facetT *facet, **facetp; + setT *facets= qh_settemp (qh num_facets); + int size, count; + + trace2((qh ferr, "qh_markkeep: only keep %d largest and/or %d most merged facets and/or min area %.2g\n", + qh KEEParea, qh KEEPmerge, qh KEEPminArea)); + FORALLfacet_(facetlist) { + if (!facet->visible && facet->good) + qh_setappend (&facets, facet); + } + size= qh_setsize (facets); + if (qh KEEParea) { + qsort (SETaddr_(facets, facetT), size, + sizeof (facetT *), qh_compare_facetarea); + if ((count= size - qh KEEParea) > 0) { + FOREACHfacet_(facets) { + facet->good= False; + if (--count == 0) + break; + } + } + } + if (qh KEEPmerge) { + qsort (SETaddr_(facets, facetT), size, + sizeof (facetT *), qh_compare_facetmerge); + if ((count= size - qh KEEPmerge) > 0) { + FOREACHfacet_(facets) { + facet->good= False; + if (--count == 0) + break; + } + } + } + if (qh KEEPminArea < REALmax/2) { + FOREACHfacet_(facets) { + if (!facet->isarea || facet->f.area < qh KEEPminArea) + facet->good= False; + } + } + qh_settempfree (&facets); + count= 0; + FORALLfacet_(facetlist) { + if (facet->good) + count++; + } + qh num_good= count; +} /* markkeep */ + + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="markvoronoi">-</a> + + qh_markvoronoi( facetlist, facets, printall, islower, numcenters ) + mark voronoi vertices for printing by site pairs + + returns: + temporary set of vertices indexed by pointid + islower set if printing lower hull (i.e., at least one facet is lower hull) + numcenters= total number of Voronoi vertices + bumps qh.printoutnum for vertex-at-infinity + clears all facet->seen and sets facet->seen2 + + if selected + facet->visitid= Voronoi vertex id + else if upper hull (or 'Qu' and lower hull) + facet->visitid= 0 + else + facet->visitid >= qh num_facets + + notes: + ignores qh.ATinfinity, if defined +*/ +setT *qh_markvoronoi (facetT *facetlist, setT *facets, boolT printall, boolT *islowerp, int *numcentersp) { + int numcenters=0; + facetT *facet, **facetp; + setT *vertices; + boolT islower= False; + + qh printoutnum++; + qh_clearcenters (qh_ASvoronoi); /* in case, qh_printvdiagram2 called by user */ + qh_vertexneighbors(); + vertices= qh_pointvertex(); + if (qh ATinfinity) + SETelem_(vertices, qh num_points-1)= NULL; + qh visit_id++; + maximize_(qh visit_id, (unsigned) qh num_facets); + FORALLfacet_(facetlist) { + if (printall || !qh_skipfacet (facet)) { + if (!facet->upperdelaunay) { + islower= True; + break; + } + } + } + FOREACHfacet_(facets) { + if (printall || !qh_skipfacet (facet)) { + if (!facet->upperdelaunay) { + islower= True; + break; + } + } + } + FORALLfacets { + if (facet->normal && (facet->upperdelaunay == islower)) + facet->visitid= 0; /* facetlist or facets may overwrite */ + else + facet->visitid= qh visit_id; + facet->seen= False; + facet->seen2= True; + } + numcenters++; /* qh_INFINITE */ + FORALLfacet_(facetlist) { + if (printall || !qh_skipfacet (facet)) + facet->visitid= numcenters++; + } + FOREACHfacet_(facets) { + if (printall || !qh_skipfacet (facet)) + facet->visitid= numcenters++; + } + *islowerp= islower; + *numcentersp= numcenters; + trace2((qh ferr, "qh_markvoronoi: islower %d numcenters %d\n", islower, numcenters)); + return vertices; +} /* markvoronoi */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="order_vertexneighbors">-</a> + + qh_order_vertexneighbors( vertex ) + order facet neighbors of a 2-d or 3-d vertex by adjacency + + notes: + does not orient the neighbors + + design: + initialize a new neighbor set with the first facet in vertex->neighbors + while vertex->neighbors non-empty + select next neighbor in the previous facet's neighbor set + set vertex->neighbors to the new neighbor set +*/ +void qh_order_vertexneighbors(vertexT *vertex) { + setT *newset; + facetT *facet, *neighbor, **neighborp; + + trace4((qh ferr, "qh_order_vertexneighbors: order neighbors of v%d for 3-d\n", vertex->id)); + newset= qh_settemp (qh_setsize (vertex->neighbors)); + facet= (facetT*)qh_setdellast (vertex->neighbors); + qh_setappend (&newset, facet); + while (qh_setsize (vertex->neighbors)) { + FOREACHneighbor_(vertex) { + if (qh_setin (facet->neighbors, neighbor)) { + qh_setdel(vertex->neighbors, neighbor); + qh_setappend (&newset, neighbor); + facet= neighbor; + break; + } + } + if (!neighbor) { + fprintf (qh ferr, "qhull internal error (qh_order_vertexneighbors): no neighbor of v%d for f%d\n", + vertex->id, facet->id); + qh_errexit (qh_ERRqhull, facet, NULL); + } + } + qh_setfree (&vertex->neighbors); + qh_settemppop (); + vertex->neighbors= newset; +} /* order_vertexneighbors */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printafacet">-</a> + + qh_printafacet( fp, format, facet, printall ) + print facet to fp in given output format (see qh.PRINTout) + + returns: + nop if !printall and qh_skipfacet() + nop if visible facet and NEWfacets and format != PRINTfacets + must match qh_countfacets + + notes + preserves qh.visit_id + facet->normal may be null if PREmerge/MERGEexact and STOPcone before merge + + see + qh_printbegin() and qh_printend() + + design: + test for printing facet + call appropriate routine for format + or output results directly +*/ +void qh_printafacet(FILE *fp, int format, facetT *facet, boolT printall) { + realT color[4], offset, dist, outerplane, innerplane; + boolT zerodiv; + coordT *point, *normp, *coordp, **pointp, *feasiblep; + int k; + vertexT *vertex, **vertexp; + facetT *neighbor, **neighborp; + + if (!printall && qh_skipfacet (facet)) + return; + if (facet->visible && qh NEWfacets && format != qh_PRINTfacets) + return; + qh printoutnum++; + switch (format) { + case qh_PRINTarea: + if (facet->isarea) { + fprintf (fp, qh_REAL_1, facet->f.area); + fprintf (fp, "\n"); + }else + fprintf (fp, "0\n"); + break; + case qh_PRINTcoplanars: + fprintf (fp, "%d", qh_setsize (facet->coplanarset)); + FOREACHpoint_(facet->coplanarset) + fprintf (fp, " %d", qh_pointid (point)); + fprintf (fp, "\n"); + break; + case qh_PRINTcentrums: + qh_printcenter (fp, format, NULL, facet); + break; + case qh_PRINTfacets: + qh_printfacet (fp, facet); + break; + case qh_PRINTfacets_xridge: + qh_printfacetheader (fp, facet); + break; + case qh_PRINTgeom: /* either 2 , 3, or 4-d by qh_printbegin */ + if (!facet->normal) + break; + for (k= qh hull_dim; k--; ) { + color[k]= (facet->normal[k]+1.0)/2.0; + maximize_(color[k], -1.0); + minimize_(color[k], +1.0); + } + qh_projectdim3 (color, color); + if (qh PRINTdim != qh hull_dim) + qh_normalize2 (color, 3, True, NULL, NULL); + if (qh hull_dim <= 2) + qh_printfacet2geom (fp, facet, color); + else if (qh hull_dim == 3) { + if (facet->simplicial) + qh_printfacet3geom_simplicial (fp, facet, color); + else + qh_printfacet3geom_nonsimplicial (fp, facet, color); + }else { + if (facet->simplicial) + qh_printfacet4geom_simplicial (fp, facet, color); + else + qh_printfacet4geom_nonsimplicial (fp, facet, color); + } + break; + case qh_PRINTids: + fprintf (fp, "%d\n", facet->id); + break; + case qh_PRINTincidences: + case qh_PRINToff: + case qh_PRINTtriangles: + if (qh hull_dim == 3 && format != qh_PRINTtriangles) + qh_printfacet3vertex (fp, facet, format); + else if (facet->simplicial || qh hull_dim == 2 || format == qh_PRINToff) + qh_printfacetNvertex_simplicial (fp, facet, format); + else + qh_printfacetNvertex_nonsimplicial (fp, facet, qh printoutvar++, format); + break; + case qh_PRINTinner: + qh_outerinner (facet, NULL, &innerplane); + offset= facet->offset - innerplane; + goto LABELprintnorm; + break; /* prevent warning */ + case qh_PRINTmerges: + fprintf (fp, "%d\n", facet->nummerge); + break; + case qh_PRINTnormals: + offset= facet->offset; + goto LABELprintnorm; + break; /* prevent warning */ + case qh_PRINTouter: + qh_outerinner (facet, &outerplane, NULL); + offset= facet->offset - outerplane; + LABELprintnorm: + if (!facet->normal) { + fprintf (fp, "no normal for facet f%d\n", facet->id); + break; + } + if (qh CDDoutput) { + fprintf (fp, qh_REAL_1, -offset); + for (k=0; k < qh hull_dim; k++) + fprintf (fp, qh_REAL_1, -facet->normal[k]); + }else { + for (k=0; k < qh hull_dim; k++) + fprintf (fp, qh_REAL_1, facet->normal[k]); + fprintf (fp, qh_REAL_1, offset); + } + fprintf (fp, "\n"); + break; + case qh_PRINTmathematica: /* either 2 or 3-d by qh_printbegin */ + if (qh hull_dim == 2) + qh_printfacet2math (fp, facet, qh printoutvar++); + else + qh_printfacet3math (fp, facet, qh printoutvar++); + break; + case qh_PRINTneighbors: + fprintf (fp, "%d", qh_setsize (facet->neighbors)); + FOREACHneighbor_(facet) + fprintf (fp, " %d", + neighbor->visitid ? neighbor->visitid - 1: - neighbor->id); + fprintf (fp, "\n"); + break; + case qh_PRINTpointintersect: + if (!qh feasible_point) { + fprintf (fp, "qhull input error (qh_printafacet): option 'Fp' needs qh feasible_point\n"); + qh_errexit( qh_ERRinput, NULL, NULL); + } + if (facet->offset > 0) + goto LABELprintinfinite; + point= coordp= (coordT*)qh_memalloc (qh normal_size); + normp= facet->normal; + feasiblep= qh feasible_point; + if (facet->offset < -qh MINdenom) { + for (k= qh hull_dim; k--; ) + *(coordp++)= (*(normp++) / - facet->offset) + *(feasiblep++); + }else { + for (k= qh hull_dim; k--; ) { + *(coordp++)= qh_divzero (*(normp++), facet->offset, qh MINdenom_1, + &zerodiv) + *(feasiblep++); + if (zerodiv) { + qh_memfree (point, qh normal_size); + goto LABELprintinfinite; + } + } + } + qh_printpoint (fp, NULL, point); + qh_memfree (point, qh normal_size); + break; + LABELprintinfinite: + for (k= qh hull_dim; k--; ) + fprintf (fp, qh_REAL_1, qh_INFINITE); + fprintf (fp, "\n"); + break; + case qh_PRINTpointnearest: + FOREACHpoint_(facet->coplanarset) { + int id, id2; + vertex= qh_nearvertex (facet, point, &dist); + id= qh_pointid (vertex->point); + id2= qh_pointid (point); + fprintf (fp, "%d %d %d " qh_REAL_1 "\n", id, id2, facet->id, dist); + } + break; + case qh_PRINTpoints: /* VORONOI only by qh_printbegin */ + if (qh CDDoutput) + fprintf (fp, "1 "); + qh_printcenter (fp, format, NULL, facet); + break; + case qh_PRINTvertices: + fprintf (fp, "%d", qh_setsize (facet->vertices)); + FOREACHvertex_(facet->vertices) + fprintf (fp, " %d", qh_pointid (vertex->point)); + fprintf (fp, "\n"); + break; + } +} /* printafacet */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printbegin">-</a> + + qh_printbegin( ) + prints header for all output formats + + returns: + checks for valid format + + notes: + uses qh.visit_id for 3/4off + changes qh.interior_point if printing centrums + qh_countfacets clears facet->visitid for non-good facets + + see + qh_printend() and qh_printafacet() + + design: + count facets and related statistics + print header for format +*/ +void qh_printbegin (FILE *fp, int format, facetT *facetlist, setT *facets, boolT printall) { + int numfacets, numsimplicial, numridges, totneighbors, numcoplanars, numtricoplanars; + int i, num; + facetT *facet, **facetp; + vertexT *vertex, **vertexp; + setT *vertices; + pointT *point, **pointp, *pointtemp; + + qh printoutnum= 0; + qh_countfacets (facetlist, facets, printall, &numfacets, &numsimplicial, + &totneighbors, &numridges, &numcoplanars, &numtricoplanars); + switch (format) { + case qh_PRINTnone: + break; + case qh_PRINTarea: + fprintf (fp, "%d\n", numfacets); + break; + case qh_PRINTcoplanars: + fprintf (fp, "%d\n", numfacets); + break; + case qh_PRINTcentrums: + if (qh CENTERtype == qh_ASnone) + qh_clearcenters (qh_AScentrum); + fprintf (fp, "%d\n%d\n", qh hull_dim, numfacets); + break; + case qh_PRINTfacets: + case qh_PRINTfacets_xridge: + if (facetlist) + qh_printvertexlist (fp, "Vertices and facets:\n", facetlist, facets, printall); + break; + case qh_PRINTgeom: + if (qh hull_dim > 4) /* qh_initqhull_globals also checks */ + goto LABELnoformat; + if (qh VORONOI && qh hull_dim > 3) /* PRINTdim == DROPdim == hull_dim-1 */ + goto LABELnoformat; + if (qh hull_dim == 2 && (qh PRINTridges || qh DOintersections)) + fprintf (qh ferr, "qhull warning: output for ridges and intersections not implemented in 2-d\n"); + if (qh hull_dim == 4 && (qh PRINTinner || qh PRINTouter || + (qh PRINTdim == 4 && qh PRINTcentrums))) + fprintf (qh ferr, "qhull warning: output for outer/inner planes and centrums not implemented in 4-d\n"); + if (qh PRINTdim == 4 && (qh PRINTspheres)) + fprintf (qh ferr, "qhull warning: output for vertices not implemented in 4-d\n"); + if (qh PRINTdim == 4 && qh DOintersections && qh PRINTnoplanes) + fprintf (qh ferr, "qhull warning: 'Gnh' generates no output in 4-d\n"); + if (qh PRINTdim == 2) { + fprintf(fp, "{appearance {linewidth 3} LIST # %s | %s\n", + qh rbox_command, qh qhull_command); + }else if (qh PRINTdim == 3) { + fprintf(fp, "{appearance {+edge -evert linewidth 2} LIST # %s | %s\n", + qh rbox_command, qh qhull_command); + }else if (qh PRINTdim == 4) { + qh visit_id++; + num= 0; + FORALLfacet_(facetlist) /* get number of ridges to be printed */ + qh_printend4geom (NULL, facet, &num, printall); + FOREACHfacet_(facets) + qh_printend4geom (NULL, facet, &num, printall); + qh ridgeoutnum= num; + qh printoutvar= 0; /* counts number of ridges in output */ + fprintf (fp, "LIST # %s | %s\n", qh rbox_command, qh qhull_command); + } + if (qh PRINTdots) { + qh printoutnum++; + num= qh num_points + qh_setsize (qh other_points); + if (qh DELAUNAY && qh ATinfinity) + num--; + if (qh PRINTdim == 4) + fprintf (fp, "4VECT %d %d 1\n", num, num); + else + fprintf (fp, "VECT %d %d 1\n", num, num); + for (i= num; i--; ) { + if (i % 20 == 0) + fprintf (fp, "\n"); + fprintf (fp, "1 "); + } + fprintf (fp, "# 1 point per line\n1 "); + for (i= num-1; i--; ) { + if (i % 20 == 0) + fprintf (fp, "\n"); + fprintf (fp, "0 "); + } + fprintf (fp, "# 1 color for all\n"); + FORALLpoints { + if (!qh DELAUNAY || !qh ATinfinity || qh_pointid(point) != qh num_points-1) { + if (qh PRINTdim == 4) + qh_printpoint (fp, NULL, point); + else + qh_printpoint3 (fp, point); + } + } + FOREACHpoint_(qh other_points) { + if (qh PRINTdim == 4) + qh_printpoint (fp, NULL, point); + else + qh_printpoint3 (fp, point); + } + fprintf (fp, "0 1 1 1 # color of points\n"); + } + if (qh PRINTdim == 4 && !qh PRINTnoplanes) + /* 4dview loads up multiple 4OFF objects slowly */ + fprintf(fp, "4OFF %d %d 1\n", 3*qh ridgeoutnum, qh ridgeoutnum); + qh PRINTcradius= 2 * qh DISTround; /* include test DISTround */ + if (qh PREmerge) { + maximize_(qh PRINTcradius, qh premerge_centrum + qh DISTround); + }else if (qh POSTmerge) + maximize_(qh PRINTcradius, qh postmerge_centrum + qh DISTround); + qh PRINTradius= qh PRINTcradius; + if (qh PRINTspheres + qh PRINTcoplanar) + maximize_(qh PRINTradius, qh MAXabs_coord * qh_MINradius); + if (qh premerge_cos < REALmax/2) { + maximize_(qh PRINTradius, (1- qh premerge_cos) * qh MAXabs_coord); + }else if (!qh PREmerge && qh POSTmerge && qh postmerge_cos < REALmax/2) { + maximize_(qh PRINTradius, (1- qh postmerge_cos) * qh MAXabs_coord); + } + maximize_(qh PRINTradius, qh MINvisible); + if (qh JOGGLEmax < REALmax/2) + qh PRINTradius += qh JOGGLEmax * sqrt (qh hull_dim); + if (qh PRINTdim != 4 && + (qh PRINTcoplanar || qh PRINTspheres || qh PRINTcentrums)) { + vertices= qh_facetvertices (facetlist, facets, printall); + if (qh PRINTspheres && qh PRINTdim <= 3) + qh_printspheres (fp, vertices, qh PRINTradius); + if (qh PRINTcoplanar || qh PRINTcentrums) { + qh firstcentrum= True; + if (qh PRINTcoplanar&& !qh PRINTspheres) { + FOREACHvertex_(vertices) + qh_printpointvect2 (fp, vertex->point, NULL, + qh interior_point, qh PRINTradius); + } + FORALLfacet_(facetlist) { + if (!printall && qh_skipfacet(facet)) + continue; + if (!facet->normal) + continue; + if (qh PRINTcentrums && qh PRINTdim <= 3) + qh_printcentrum (fp, facet, qh PRINTcradius); + if (!qh PRINTcoplanar) + continue; + FOREACHpoint_(facet->coplanarset) + qh_printpointvect2 (fp, point, facet->normal, NULL, qh PRINTradius); + FOREACHpoint_(facet->outsideset) + qh_printpointvect2 (fp, point, facet->normal, NULL, qh PRINTradius); + } + FOREACHfacet_(facets) { + if (!printall && qh_skipfacet(facet)) + continue; + if (!facet->normal) + continue; + if (qh PRINTcentrums && qh PRINTdim <= 3) + qh_printcentrum (fp, facet, qh PRINTcradius); + if (!qh PRINTcoplanar) + continue; + FOREACHpoint_(facet->coplanarset) + qh_printpointvect2 (fp, point, facet->normal, NULL, qh PRINTradius); + FOREACHpoint_(facet->outsideset) + qh_printpointvect2 (fp, point, facet->normal, NULL, qh PRINTradius); + } + } + qh_settempfree (&vertices); + } + qh visit_id++; /* for printing hyperplane intersections */ + break; + case qh_PRINTids: + fprintf (fp, "%d\n", numfacets); + break; + case qh_PRINTincidences: + if (qh VORONOI && qh PRINTprecision) + fprintf (qh ferr, "qhull warning: writing Delaunay. Use 'p' or 'o' for Voronoi centers\n"); + qh printoutvar= qh vertex_id; /* centrum id for non-simplicial facets */ + if (qh hull_dim <= 3) + fprintf(fp, "%d\n", numfacets); + else + fprintf(fp, "%d\n", numsimplicial+numridges); + break; + case qh_PRINTinner: + case qh_PRINTnormals: + case qh_PRINTouter: + if (qh CDDoutput) + fprintf (fp, "%s | %s\nbegin\n %d %d real\n", qh rbox_command, + qh qhull_command, numfacets, qh hull_dim+1); + else + fprintf (fp, "%d\n%d\n", qh hull_dim+1, numfacets); + break; + case qh_PRINTmathematica: + if (qh hull_dim > 3) /* qh_initbuffers also checks */ + goto LABELnoformat; + if (qh VORONOI) + fprintf (qh ferr, "qhull warning: output is the Delaunay triangulation\n"); + fprintf(fp, "{\n"); + qh printoutvar= 0; /* counts number of facets for notfirst */ + break; + case qh_PRINTmerges: + fprintf (fp, "%d\n", numfacets); + break; + case qh_PRINTpointintersect: + fprintf (fp, "%d\n%d\n", qh hull_dim, numfacets); + break; + case qh_PRINTneighbors: + fprintf (fp, "%d\n", numfacets); + break; + case qh_PRINToff: + case qh_PRINTtriangles: + if (qh VORONOI) + goto LABELnoformat; + num = qh hull_dim; + if (format == qh_PRINToff || qh hull_dim == 2) + fprintf (fp, "%d\n%d %d %d\n", num, + qh num_points+qh_setsize (qh other_points), numfacets, totneighbors/2); + else { /* qh_PRINTtriangles */ + qh printoutvar= qh num_points+qh_setsize (qh other_points); /* first centrum */ + if (qh DELAUNAY) + num--; /* drop last dimension */ + fprintf (fp, "%d\n%d %d %d\n", num, qh printoutvar + + numfacets - numsimplicial, numsimplicial + numridges, totneighbors/2); + } + FORALLpoints + qh_printpointid (qh fout, NULL, num, point, -1); + FOREACHpoint_(qh other_points) + qh_printpointid (qh fout, NULL, num, point, -1); + if (format == qh_PRINTtriangles && qh hull_dim > 2) { + FORALLfacets { + if (!facet->simplicial && facet->visitid) + qh_printcenter (qh fout, format, NULL, facet); + } + } + break; + case qh_PRINTpointnearest: + fprintf (fp, "%d\n", numcoplanars); + break; + case qh_PRINTpoints: + if (!qh VORONOI) + goto LABELnoformat; + if (qh CDDoutput) + fprintf (fp, "%s | %s\nbegin\n%d %d real\n", qh rbox_command, + qh qhull_command, numfacets, qh hull_dim); + else + fprintf (fp, "%d\n%d\n", qh hull_dim-1, numfacets); + break; + case qh_PRINTvertices: + fprintf (fp, "%d\n", numfacets); + break; + case qh_PRINTsummary: + default: + LABELnoformat: + fprintf (qh ferr, "qhull internal error (qh_printbegin): can not use this format for dimension %d\n", + qh hull_dim); + qh_errexit (qh_ERRqhull, NULL, NULL); + } +} /* printbegin */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printcenter">-</a> + + qh_printcenter( fp, string, facet ) + print facet->center as centrum or Voronoi center + string may be NULL. Don't include '%' codes. + nop if qh CENTERtype neither CENTERvoronoi nor CENTERcentrum + if upper envelope of Delaunay triangulation and point at-infinity + prints qh_INFINITE instead; + + notes: + defines facet->center if needed + if format=PRINTgeom, adds a 0 if would otherwise be 2-d +*/ +void qh_printcenter (FILE *fp, int format, char *string, facetT *facet) { + int k, num; + + if (qh CENTERtype != qh_ASvoronoi && qh CENTERtype != qh_AScentrum) + return; + if (string) + fprintf (fp, string, facet->id); + if (qh CENTERtype == qh_ASvoronoi) { + num= qh hull_dim-1; + if (!facet->normal || !facet->upperdelaunay || !qh ATinfinity) { + if (!facet->center) + facet->center= qh_facetcenter (facet->vertices); + for (k=0; k < num; k++) + fprintf (fp, qh_REAL_1, facet->center[k]); + }else { + for (k=0; k < num; k++) + fprintf (fp, qh_REAL_1, qh_INFINITE); + } + }else /* qh CENTERtype == qh_AScentrum */ { + num= qh hull_dim; + if (format == qh_PRINTtriangles && qh DELAUNAY) + num--; + if (!facet->center) + facet->center= qh_getcentrum (facet); + for (k=0; k < num; k++) + fprintf (fp, qh_REAL_1, facet->center[k]); + } + if (format == qh_PRINTgeom && num == 2) + fprintf (fp, " 0\n"); + else + fprintf (fp, "\n"); +} /* printcenter */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printcentrum">-</a> + + qh_printcentrum( fp, facet, radius ) + print centrum for a facet in OOGL format + radius defines size of centrum + 2-d or 3-d only + + returns: + defines facet->center if needed +*/ +void qh_printcentrum (FILE *fp, facetT *facet, realT radius) { + pointT *centrum, *projpt; + boolT tempcentrum= False; + realT xaxis[4], yaxis[4], normal[4], dist; + realT green[3]={0, 1, 0}; + vertexT *apex; + int k; + + if (qh CENTERtype == qh_AScentrum) { + if (!facet->center) + facet->center= qh_getcentrum (facet); + centrum= facet->center; + }else { + centrum= qh_getcentrum (facet); + tempcentrum= True; + } + fprintf (fp, "{appearance {-normal -edge normscale 0} "); + if (qh firstcentrum) { + qh firstcentrum= False; + fprintf (fp, "{INST geom { define centrum CQUAD # f%d\n\ +-0.3 -0.3 0.0001 0 0 1 1\n\ + 0.3 -0.3 0.0001 0 0 1 1\n\ + 0.3 0.3 0.0001 0 0 1 1\n\ +-0.3 0.3 0.0001 0 0 1 1 } transform { \n", facet->id); + }else + fprintf (fp, "{INST geom { : centrum } transform { # f%d\n", facet->id); + apex= SETfirstt_(facet->vertices, vertexT); + qh_distplane(apex->point, facet, &dist); + projpt= qh_projectpoint(apex->point, facet, dist); + for (k= qh hull_dim; k--; ) { + xaxis[k]= projpt[k] - centrum[k]; + normal[k]= facet->normal[k]; + } + if (qh hull_dim == 2) { + xaxis[2]= 0; + normal[2]= 0; + }else if (qh hull_dim == 4) { + qh_projectdim3 (xaxis, xaxis); + qh_projectdim3 (normal, normal); + qh_normalize2 (normal, qh PRINTdim, True, NULL, NULL); + } + qh_crossproduct (3, xaxis, normal, yaxis); + fprintf (fp, "%8.4g %8.4g %8.4g 0\n", xaxis[0], xaxis[1], xaxis[2]); + fprintf (fp, "%8.4g %8.4g %8.4g 0\n", yaxis[0], yaxis[1], yaxis[2]); + fprintf (fp, "%8.4g %8.4g %8.4g 0\n", normal[0], normal[1], normal[2]); + qh_printpoint3 (fp, centrum); + fprintf (fp, "1 }}}\n"); + qh_memfree (projpt, qh normal_size); + qh_printpointvect (fp, centrum, facet->normal, NULL, radius, green); + if (tempcentrum) + qh_memfree (centrum, qh normal_size); +} /* printcentrum */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printend">-</a> + + qh_printend( fp, format ) + prints trailer for all output formats + + see: + qh_printbegin() and qh_printafacet() + +*/ +void qh_printend (FILE *fp, int format, facetT *facetlist, setT *facets, boolT printall) { + int num; + facetT *facet, **facetp; + + if (!qh printoutnum) + fprintf (qh ferr, "qhull warning: no facets printed\n"); + switch (format) { + case qh_PRINTgeom: + if (qh hull_dim == 4 && qh DROPdim < 0 && !qh PRINTnoplanes) { + qh visit_id++; + num= 0; + FORALLfacet_(facetlist) + qh_printend4geom (fp, facet,&num, printall); + FOREACHfacet_(facets) + qh_printend4geom (fp, facet, &num, printall); + if (num != qh ridgeoutnum || qh printoutvar != qh ridgeoutnum) { + fprintf (qh ferr, "qhull internal error (qh_printend): number of ridges %d != number printed %d and at end %d\n", qh ridgeoutnum, qh printoutvar, num); + qh_errexit (qh_ERRqhull, NULL, NULL); + } + }else + fprintf(fp, "}\n"); + break; + case qh_PRINTinner: + case qh_PRINTnormals: + case qh_PRINTouter: + if (qh CDDoutput) + fprintf (fp, "end\n"); + break; + case qh_PRINTmathematica: + fprintf(fp, "}\n"); + break; + case qh_PRINTpoints: + if (qh CDDoutput) + fprintf (fp, "end\n"); + break; + } +} /* printend */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printend4geom">-</a> + + qh_printend4geom( fp, facet, numridges, printall ) + helper function for qh_printbegin/printend + + returns: + number of printed ridges + + notes: + just counts printed ridges if fp=NULL + uses facet->visitid + must agree with qh_printfacet4geom... + + design: + computes color for facet from its normal + prints each ridge of facet +*/ +void qh_printend4geom (FILE *fp, facetT *facet, int *nump, boolT printall) { + realT color[3]; + int i, num= *nump; + facetT *neighbor, **neighborp; + ridgeT *ridge, **ridgep; + + if (!printall && qh_skipfacet(facet)) + return; + if (qh PRINTnoplanes || (facet->visible && qh NEWfacets)) + return; + if (!facet->normal) + return; + if (fp) { + for (i=0; i < 3; i++) { + color[i]= (facet->normal[i]+1.0)/2.0; + maximize_(color[i], -1.0); + minimize_(color[i], +1.0); + } + } + facet->visitid= qh visit_id; + if (facet->simplicial) { + FOREACHneighbor_(facet) { + if (neighbor->visitid != qh visit_id) { + if (fp) + fprintf (fp, "3 %d %d %d %8.4g %8.4g %8.4g 1 # f%d f%d\n", + 3*num, 3*num+1, 3*num+2, color[0], color[1], color[2], + facet->id, neighbor->id); + num++; + } + } + }else { + FOREACHridge_(facet->ridges) { + neighbor= otherfacet_(ridge, facet); + if (neighbor->visitid != qh visit_id) { + if (fp) + fprintf (fp, "3 %d %d %d %8.4g %8.4g %8.4g 1 #r%d f%d f%d\n", + 3*num, 3*num+1, 3*num+2, color[0], color[1], color[2], + ridge->id, facet->id, neighbor->id); + num++; + } + } + } + *nump= num; +} /* printend4geom */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printextremes">-</a> + + qh_printextremes( fp, facetlist, facets, printall ) + print extreme points for convex hulls or halfspace intersections + + notes: + #points, followed by ids, one per line + + sorted by id + same order as qh_printpoints_out if no coplanar/interior points +*/ +void qh_printextremes (FILE *fp, facetT *facetlist, setT *facets, int printall) { + setT *vertices, *points; + pointT *point; + vertexT *vertex, **vertexp; + int id; + int numpoints=0, point_i, point_n; + int allpoints= qh num_points + qh_setsize (qh other_points); + + points= qh_settemp (allpoints); + qh_setzero (points, 0, allpoints); + vertices= qh_facetvertices (facetlist, facets, printall); + FOREACHvertex_(vertices) { + id= qh_pointid (vertex->point); + if (id >= 0) { + SETelem_(points, id)= vertex->point; + numpoints++; + } + } + qh_settempfree (&vertices); + fprintf (fp, "%d\n", numpoints); + FOREACHpoint_i_(points) { + if (point) + fprintf (fp, "%d\n", point_i); + } + qh_settempfree (&points); +} /* printextremes */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printextremes_2d">-</a> + + qh_printextremes_2d( fp, facetlist, facets, printall ) + prints point ids for facets in qh_ORIENTclock order + + notes: + #points, followed by ids, one per line + if facetlist/facets are disjoint than the output includes skips + errors if facets form a loop + does not print coplanar points +*/ +void qh_printextremes_2d (FILE *fp, facetT *facetlist, setT *facets, int printall) { + int numfacets, numridges, totneighbors, numcoplanars, numsimplicial, numtricoplanars; + setT *vertices; + facetT *facet, *startfacet, *nextfacet; + vertexT *vertexA, *vertexB; + + qh_countfacets (facetlist, facets, printall, &numfacets, &numsimplicial, + &totneighbors, &numridges, &numcoplanars, &numtricoplanars); /* marks qh visit_id */ + vertices= qh_facetvertices (facetlist, facets, printall); + fprintf(fp, "%d\n", qh_setsize (vertices)); + qh_settempfree (&vertices); + if (!numfacets) + return; + facet= startfacet= facetlist ? facetlist : SETfirstt_(facets, facetT); + qh vertex_visit++; + qh visit_id++; + do { + if (facet->toporient ^ qh_ORIENTclock) { + vertexA= SETfirstt_(facet->vertices, vertexT); + vertexB= SETsecondt_(facet->vertices, vertexT); + nextfacet= SETfirstt_(facet->neighbors, facetT); + }else { + vertexA= SETsecondt_(facet->vertices, vertexT); + vertexB= SETfirstt_(facet->vertices, vertexT); + nextfacet= SETsecondt_(facet->neighbors, facetT); + } + if (facet->visitid == qh visit_id) { + fprintf(qh ferr, "qh_printextremes_2d: loop in facet list. facet %d nextfacet %d\n", + facet->id, nextfacet->id); + qh_errexit2 (qh_ERRqhull, facet, nextfacet); + } + if (facet->visitid) { + if (vertexA->visitid != qh vertex_visit) { + vertexA->visitid= qh vertex_visit; + fprintf(fp, "%d\n", qh_pointid (vertexA->point)); + } + if (vertexB->visitid != qh vertex_visit) { + vertexB->visitid= qh vertex_visit; + fprintf(fp, "%d\n", qh_pointid (vertexB->point)); + } + } + facet->visitid= qh visit_id; + facet= nextfacet; + }while (facet && facet != startfacet); +} /* printextremes_2d */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printextremes_d">-</a> + + qh_printextremes_d( fp, facetlist, facets, printall ) + print extreme points of input sites for Delaunay triangulations + + notes: + #points, followed by ids, one per line + + unordered +*/ +void qh_printextremes_d (FILE *fp, facetT *facetlist, setT *facets, int printall) { + setT *vertices; + vertexT *vertex, **vertexp; + boolT upperseen, lowerseen; + facetT *neighbor, **neighborp; + int numpoints=0; + + vertices= qh_facetvertices (facetlist, facets, printall); + qh_vertexneighbors(); + FOREACHvertex_(vertices) { + upperseen= lowerseen= False; + FOREACHneighbor_(vertex) { + if (neighbor->upperdelaunay) + upperseen= True; + else + lowerseen= True; + } + if (upperseen && lowerseen) { + vertex->seen= True; + numpoints++; + }else + vertex->seen= False; + } + fprintf (fp, "%d\n", numpoints); + FOREACHvertex_(vertices) { + if (vertex->seen) + fprintf (fp, "%d\n", qh_pointid (vertex->point)); + } + qh_settempfree (&vertices); +} /* printextremes_d */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printfacet">-</a> + + qh_printfacet( fp, facet ) + prints all fields of a facet to fp + + notes: + ridges printed in neighbor order +*/ +void qh_printfacet(FILE *fp, facetT *facet) { + + qh_printfacetheader (fp, facet); + if (facet->ridges) + qh_printfacetridges (fp, facet); +} /* printfacet */ + + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printfacet2geom">-</a> + + qh_printfacet2geom( fp, facet, color ) + print facet as part of a 2-d VECT for Geomview + + notes: + assume precise calculations in io.c with roundoff covered by qh_GEOMepsilon + mindist is calculated within io.c. maxoutside is calculated elsewhere + so a DISTround error may have occured. +*/ +void qh_printfacet2geom(FILE *fp, facetT *facet, realT color[3]) { + pointT *point0, *point1; + realT mindist, innerplane, outerplane; + int k; + + qh_facet2point (facet, &point0, &point1, &mindist); + qh_geomplanes (facet, &outerplane, &innerplane); + if (qh PRINTouter || (!qh PRINTnoplanes && !qh PRINTinner)) + qh_printfacet2geom_points(fp, point0, point1, facet, outerplane, color); + if (qh PRINTinner || (!qh PRINTnoplanes && !qh PRINTouter && + outerplane - innerplane > 2 * qh MAXabs_coord * qh_GEOMepsilon)) { + for(k= 3; k--; ) + color[k]= 1.0 - color[k]; + qh_printfacet2geom_points(fp, point0, point1, facet, innerplane, color); + } + qh_memfree (point1, qh normal_size); + qh_memfree (point0, qh normal_size); +} /* printfacet2geom */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printfacet2geom_points">-</a> + + qh_printfacet2geom_points( fp, point1, point2, facet, offset, color ) + prints a 2-d facet as a VECT with 2 points at some offset. + The points are on the facet's plane. +*/ +void qh_printfacet2geom_points(FILE *fp, pointT *point1, pointT *point2, + facetT *facet, realT offset, realT color[3]) { + pointT *p1= point1, *p2= point2; + + fprintf(fp, "VECT 1 2 1 2 1 # f%d\n", facet->id); + if (offset != 0.0) { + p1= qh_projectpoint (p1, facet, -offset); + p2= qh_projectpoint (p2, facet, -offset); + } + fprintf(fp, "%8.4g %8.4g %8.4g\n%8.4g %8.4g %8.4g\n", + p1[0], p1[1], 0.0, p2[0], p2[1], 0.0); + if (offset != 0.0) { + qh_memfree (p1, qh normal_size); + qh_memfree (p2, qh normal_size); + } + fprintf(fp, "%8.4g %8.4g %8.4g 1.0\n", color[0], color[1], color[2]); +} /* printfacet2geom_points */ + + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printfacet2math">-</a> + + qh_printfacet2math( fp, facet, notfirst ) + print 2-d Mathematica output for a facet + may be non-simplicial + + notes: + use %16.8f since Mathematica 2.2 does not handle exponential format +*/ +void qh_printfacet2math(FILE *fp, facetT *facet, int notfirst) { + pointT *point0, *point1; + realT mindist; + + qh_facet2point (facet, &point0, &point1, &mindist); + if (notfirst) + fprintf(fp, ","); + fprintf(fp, "Line[{{%16.8f, %16.8f}, {%16.8f, %16.8f}}]\n", + point0[0], point0[1], point1[0], point1[1]); + qh_memfree (point1, qh normal_size); + qh_memfree (point0, qh normal_size); +} /* printfacet2math */ + + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printfacet3geom_nonsimplicial">-</a> + + qh_printfacet3geom_nonsimplicial( fp, facet, color ) + print Geomview OFF for a 3-d nonsimplicial facet. + if DOintersections, prints ridges to unvisited neighbors (qh visit_id) + + notes + uses facet->visitid for intersections and ridges +*/ +void qh_printfacet3geom_nonsimplicial(FILE *fp, facetT *facet, realT color[3]) { + ridgeT *ridge, **ridgep; + setT *projectedpoints, *vertices; + vertexT *vertex, **vertexp, *vertexA, *vertexB; + pointT *projpt, *point, **pointp; + facetT *neighbor; + realT dist, outerplane, innerplane; + int cntvertices, k; + realT black[3]={0, 0, 0}, green[3]={0, 1, 0}; + + qh_geomplanes (facet, &outerplane, &innerplane); + vertices= qh_facet3vertex (facet); /* oriented */ + cntvertices= qh_setsize(vertices); + projectedpoints= qh_settemp(cntvertices); + FOREACHvertex_(vertices) { + zinc_(Zdistio); + qh_distplane(vertex->point, facet, &dist); + projpt= qh_projectpoint(vertex->point, facet, dist); + qh_setappend (&projectedpoints, projpt); + } + if (qh PRINTouter || (!qh PRINTnoplanes && !qh PRINTinner)) + qh_printfacet3geom_points(fp, projectedpoints, facet, outerplane, color); + if (qh PRINTinner || (!qh PRINTnoplanes && !qh PRINTouter && + outerplane - innerplane > 2 * qh MAXabs_coord * qh_GEOMepsilon)) { + for (k=3; k--; ) + color[k]= 1.0 - color[k]; + qh_printfacet3geom_points(fp, projectedpoints, facet, innerplane, color); + } + FOREACHpoint_(projectedpoints) + qh_memfree (point, qh normal_size); + qh_settempfree(&projectedpoints); + qh_settempfree(&vertices); + if ((qh DOintersections || qh PRINTridges) + && (!facet->visible || !qh NEWfacets)) { + facet->visitid= qh visit_id; + FOREACHridge_(facet->ridges) { + neighbor= otherfacet_(ridge, facet); + if (neighbor->visitid != qh visit_id) { + if (qh DOintersections) + qh_printhyperplaneintersection(fp, facet, neighbor, ridge->vertices, black); + if (qh PRINTridges) { + vertexA= SETfirstt_(ridge->vertices, vertexT); + vertexB= SETsecondt_(ridge->vertices, vertexT); + qh_printline3geom (fp, vertexA->point, vertexB->point, green); + } + } + } + } +} /* printfacet3geom_nonsimplicial */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printfacet3geom_points">-</a> + + qh_printfacet3geom_points( fp, points, facet, offset ) + prints a 3-d facet as OFF Geomview object. + offset is relative to the facet's hyperplane + Facet is determined as a list of points +*/ +void qh_printfacet3geom_points(FILE *fp, setT *points, facetT *facet, realT offset, realT color[3]) { + int k, n= qh_setsize(points), i; + pointT *point, **pointp; + setT *printpoints; + + fprintf(fp, "{ OFF %d 1 1 # f%d\n", n, facet->id); + if (offset != 0.0) { + printpoints= qh_settemp (n); + FOREACHpoint_(points) + qh_setappend (&printpoints, qh_projectpoint(point, facet, -offset)); + }else + printpoints= points; + FOREACHpoint_(printpoints) { + for (k=0; k < qh hull_dim; k++) { + if (k == qh DROPdim) + fprintf(fp, "0 "); + else + fprintf(fp, "%8.4g ", point[k]); + } + if (printpoints != points) + qh_memfree (point, qh normal_size); + fprintf (fp, "\n"); + } + if (printpoints != points) + qh_settempfree (&printpoints); + fprintf(fp, "%d ", n); + for(i= 0; i < n; i++) + fprintf(fp, "%d ", i); + fprintf(fp, "%8.4g %8.4g %8.4g 1.0 }\n", color[0], color[1], color[2]); +} /* printfacet3geom_points */ + + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printfacet3geom_simplicial">-</a> + + qh_printfacet3geom_simplicial( ) + print Geomview OFF for a 3-d simplicial facet. + + notes: + may flip color + uses facet->visitid for intersections and ridges + + assume precise calculations in io.c with roundoff covered by qh_GEOMepsilon + innerplane may be off by qh DISTround. Maxoutside is calculated elsewhere + so a DISTround error may have occured. +*/ +void qh_printfacet3geom_simplicial(FILE *fp, facetT *facet, realT color[3]) { + setT *points, *vertices; + vertexT *vertex, **vertexp, *vertexA, *vertexB; + facetT *neighbor, **neighborp; + realT outerplane, innerplane; + realT black[3]={0, 0, 0}, green[3]={0, 1, 0}; + int k; + + qh_geomplanes (facet, &outerplane, &innerplane); + vertices= qh_facet3vertex (facet); + points= qh_settemp (qh TEMPsize); + FOREACHvertex_(vertices) + qh_setappend(&points, vertex->point); + if (qh PRINTouter || (!qh PRINTnoplanes && !qh PRINTinner)) + qh_printfacet3geom_points(fp, points, facet, outerplane, color); + if (qh PRINTinner || (!qh PRINTnoplanes && !qh PRINTouter && + outerplane - innerplane > 2 * qh MAXabs_coord * qh_GEOMepsilon)) { + for (k= 3; k--; ) + color[k]= 1.0 - color[k]; + qh_printfacet3geom_points(fp, points, facet, innerplane, color); + } + qh_settempfree(&points); + qh_settempfree(&vertices); + if ((qh DOintersections || qh PRINTridges) + && (!facet->visible || !qh NEWfacets)) { + facet->visitid= qh visit_id; + FOREACHneighbor_(facet) { + if (neighbor->visitid != qh visit_id) { + vertices= qh_setnew_delnthsorted (facet->vertices, qh hull_dim, + SETindex_(facet->neighbors, neighbor), 0); + if (qh DOintersections) + qh_printhyperplaneintersection(fp, facet, neighbor, vertices, black); + if (qh PRINTridges) { + vertexA= SETfirstt_(vertices, vertexT); + vertexB= SETsecondt_(vertices, vertexT); + qh_printline3geom (fp, vertexA->point, vertexB->point, green); + } + qh_setfree(&vertices); + } + } + } +} /* printfacet3geom_simplicial */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printfacet3math">-</a> + + qh_printfacet3math( fp, facet, notfirst ) + print 3-d Mathematica output for a facet + + notes: + may be non-simplicial + use %16.8f since Mathematica 2.2 does not handle exponential format +*/ +void qh_printfacet3math (FILE *fp, facetT *facet, int notfirst) { + vertexT *vertex, **vertexp; + setT *points, *vertices; + pointT *point, **pointp; + boolT firstpoint= True; + realT dist; + + if (notfirst) + fprintf(fp, ",\n"); + vertices= qh_facet3vertex (facet); + points= qh_settemp (qh_setsize (vertices)); + FOREACHvertex_(vertices) { + zinc_(Zdistio); + qh_distplane(vertex->point, facet, &dist); + point= qh_projectpoint(vertex->point, facet, dist); + qh_setappend (&points, point); + } + fprintf(fp, "Polygon[{"); + FOREACHpoint_(points) { + if (firstpoint) + firstpoint= False; + else + fprintf(fp, ",\n"); + fprintf(fp, "{%16.8f, %16.8f, %16.8f}", point[0], point[1], point[2]); + } + FOREACHpoint_(points) + qh_memfree (point, qh normal_size); + qh_settempfree(&points); + qh_settempfree(&vertices); + fprintf(fp, "}]"); +} /* printfacet3math */ + + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printfacet3vertex">-</a> + + qh_printfacet3vertex( fp, facet, format ) + print vertices in a 3-d facet as point ids + + notes: + prints number of vertices first if format == qh_PRINToff + the facet may be non-simplicial +*/ +void qh_printfacet3vertex(FILE *fp, facetT *facet, int format) { + vertexT *vertex, **vertexp; + setT *vertices; + + vertices= qh_facet3vertex (facet); + if (format == qh_PRINToff) + fprintf (fp, "%d ", qh_setsize (vertices)); + FOREACHvertex_(vertices) + fprintf (fp, "%d ", qh_pointid(vertex->point)); + fprintf (fp, "\n"); + qh_settempfree(&vertices); +} /* printfacet3vertex */ + + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printfacet4geom_nonsimplicial">-</a> + + qh_printfacet4geom_nonsimplicial( ) + print Geomview 4OFF file for a 4d nonsimplicial facet + prints all ridges to unvisited neighbors (qh.visit_id) + if qh.DROPdim + prints in OFF format + + notes: + must agree with printend4geom() +*/ +void qh_printfacet4geom_nonsimplicial(FILE *fp, facetT *facet, realT color[3]) { + facetT *neighbor; + ridgeT *ridge, **ridgep; + vertexT *vertex, **vertexp; + pointT *point; + int k; + realT dist; + + facet->visitid= qh visit_id; + if (qh PRINTnoplanes || (facet->visible && qh NEWfacets)) + return; + FOREACHridge_(facet->ridges) { + neighbor= otherfacet_(ridge, facet); + if (neighbor->visitid == qh visit_id) + continue; + if (qh PRINTtransparent && !neighbor->good) + continue; + if (qh DOintersections) + qh_printhyperplaneintersection(fp, facet, neighbor, ridge->vertices, color); + else { + if (qh DROPdim >= 0) + fprintf(fp, "OFF 3 1 1 # f%d\n", facet->id); + else { + qh printoutvar++; + fprintf (fp, "# r%d between f%d f%d\n", ridge->id, facet->id, neighbor->id); + } + FOREACHvertex_(ridge->vertices) { + zinc_(Zdistio); + qh_distplane(vertex->point,facet, &dist); + point=qh_projectpoint(vertex->point,facet, dist); + for(k= 0; k < qh hull_dim; k++) { + if (k != qh DROPdim) + fprintf(fp, "%8.4g ", point[k]); + } + fprintf (fp, "\n"); + qh_memfree (point, qh normal_size); + } + if (qh DROPdim >= 0) + fprintf(fp, "3 0 1 2 %8.4g %8.4g %8.4g\n", color[0], color[1], color[2]); + } + } +} /* printfacet4geom_nonsimplicial */ + + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printfacet4geom_simplicial">-</a> + + qh_printfacet4geom_simplicial( fp, facet, color ) + print Geomview 4OFF file for a 4d simplicial facet + prints triangles for unvisited neighbors (qh.visit_id) + + notes: + must agree with printend4geom() +*/ +void qh_printfacet4geom_simplicial(FILE *fp, facetT *facet, realT color[3]) { + setT *vertices; + facetT *neighbor, **neighborp; + vertexT *vertex, **vertexp; + int k; + + facet->visitid= qh visit_id; + if (qh PRINTnoplanes || (facet->visible && qh NEWfacets)) + return; + FOREACHneighbor_(facet) { + if (neighbor->visitid == qh visit_id) + continue; + if (qh PRINTtransparent && !neighbor->good) + continue; + vertices= qh_setnew_delnthsorted (facet->vertices, qh hull_dim, + SETindex_(facet->neighbors, neighbor), 0); + if (qh DOintersections) + qh_printhyperplaneintersection(fp, facet, neighbor, vertices, color); + else { + if (qh DROPdim >= 0) + fprintf(fp, "OFF 3 1 1 # ridge between f%d f%d\n", + facet->id, neighbor->id); + else { + qh printoutvar++; + fprintf (fp, "# ridge between f%d f%d\n", facet->id, neighbor->id); + } + FOREACHvertex_(vertices) { + for(k= 0; k < qh hull_dim; k++) { + if (k != qh DROPdim) + fprintf(fp, "%8.4g ", vertex->point[k]); + } + fprintf (fp, "\n"); + } + if (qh DROPdim >= 0) + fprintf(fp, "3 0 1 2 %8.4g %8.4g %8.4g\n", color[0], color[1], color[2]); + } + qh_setfree(&vertices); + } +} /* printfacet4geom_simplicial */ + + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printfacetNvertex_nonsimplicial">-</a> + + qh_printfacetNvertex_nonsimplicial( fp, facet, id, format ) + print vertices for an N-d non-simplicial facet + triangulates each ridge to the id +*/ +void qh_printfacetNvertex_nonsimplicial(FILE *fp, facetT *facet, int id, int format) { + vertexT *vertex, **vertexp; + ridgeT *ridge, **ridgep; + + if (facet->visible && qh NEWfacets) + return; + FOREACHridge_(facet->ridges) { + if (format == qh_PRINTtriangles) + fprintf(fp, "%d ", qh hull_dim); + fprintf(fp, "%d ", id); + if ((ridge->top == facet) ^ qh_ORIENTclock) { + FOREACHvertex_(ridge->vertices) + fprintf(fp, "%d ", qh_pointid(vertex->point)); + }else { + FOREACHvertexreverse12_(ridge->vertices) + fprintf(fp, "%d ", qh_pointid(vertex->point)); + } + fprintf(fp, "\n"); + } +} /* printfacetNvertex_nonsimplicial */ + + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printfacetNvertex_simplicial">-</a> + + qh_printfacetNvertex_simplicial( fp, facet, format ) + print vertices for an N-d simplicial facet + prints vertices for non-simplicial facets + 2-d facets (orientation preserved by qh_mergefacet2d) + PRINToff ('o') for 4-d and higher +*/ +void qh_printfacetNvertex_simplicial(FILE *fp, facetT *facet, int format) { + vertexT *vertex, **vertexp; + + if (format == qh_PRINToff || format == qh_PRINTtriangles) + fprintf (fp, "%d ", qh_setsize (facet->vertices)); + if ((facet->toporient ^ qh_ORIENTclock) + || (qh hull_dim > 2 && !facet->simplicial)) { + FOREACHvertex_(facet->vertices) + fprintf(fp, "%d ", qh_pointid(vertex->point)); + }else { + FOREACHvertexreverse12_(facet->vertices) + fprintf(fp, "%d ", qh_pointid(vertex->point)); + } + fprintf(fp, "\n"); +} /* printfacetNvertex_simplicial */ + + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printfacetheader">-</a> + + qh_printfacetheader( fp, facet ) + prints header fields of a facet to fp + + notes: + for 'f' output and debugging +*/ +void qh_printfacetheader(FILE *fp, facetT *facet) { + pointT *point, **pointp, *furthest; + facetT *neighbor, **neighborp; + realT dist; + + if (facet == qh_MERGEridge) { + fprintf (fp, " MERGEridge\n"); + return; + }else if (facet == qh_DUPLICATEridge) { + fprintf (fp, " DUPLICATEridge\n"); + return; + }else if (!facet) { + fprintf (fp, " NULLfacet\n"); + return; + } + qh old_randomdist= qh RANDOMdist; + qh RANDOMdist= False; + fprintf(fp, "- f%d\n", facet->id); + fprintf(fp, " - flags:"); + if (facet->toporient) + fprintf(fp, " top"); + else + fprintf(fp, " bottom"); + if (facet->simplicial) + fprintf(fp, " simplicial"); + if (facet->tricoplanar) + fprintf(fp, " tricoplanar"); + if (facet->upperdelaunay) + fprintf(fp, " upperDelaunay"); + if (facet->visible) + fprintf(fp, " visible"); + if (facet->newfacet) + fprintf(fp, " new"); + if (facet->tested) + fprintf(fp, " tested"); + if (!facet->good) + fprintf(fp, " notG"); + if (facet->seen) + fprintf(fp, " seen"); + if (facet->coplanar) + fprintf(fp, " coplanar"); + if (facet->mergehorizon) + fprintf(fp, " mergehorizon"); + if (facet->keepcentrum) + fprintf(fp, " keepcentrum"); + if (facet->dupridge) + fprintf(fp, " dupridge"); + if (facet->mergeridge && !facet->mergeridge2) + fprintf(fp, " mergeridge1"); + if (facet->mergeridge2) + fprintf(fp, " mergeridge2"); + if (facet->newmerge) + fprintf(fp, " newmerge"); + if (facet->flipped) + fprintf(fp, " flipped"); + if (facet->notfurthest) + fprintf(fp, " notfurthest"); + if (facet->degenerate) + fprintf(fp, " degenerate"); + if (facet->redundant) + fprintf(fp, " redundant"); + fprintf(fp, "\n"); + if (facet->isarea) + fprintf(fp, " - area: %2.2g\n", facet->f.area); + else if (qh NEWfacets && facet->visible && facet->f.replace) + fprintf(fp, " - replacement: f%d\n", facet->f.replace->id); + else if (facet->newfacet) { + if (facet->f.samecycle && facet->f.samecycle != facet) + fprintf(fp, " - shares same visible/horizon as f%d\n", facet->f.samecycle->id); + }else if (facet->tricoplanar /* !isarea */) { + if (facet->f.triowner) + fprintf(fp, " - owner of normal & centrum is facet f%d\n", facet->f.triowner->id); + }else if (facet->f.newcycle) + fprintf(fp, " - was horizon to f%d\n", facet->f.newcycle->id); + if (facet->nummerge) + fprintf(fp, " - merges: %d\n", facet->nummerge); + qh_printpointid(fp, " - normal: ", qh hull_dim, facet->normal, -1); + fprintf(fp, " - offset: %10.7g\n", facet->offset); + if (qh CENTERtype == qh_ASvoronoi || facet->center) + qh_printcenter (fp, qh_PRINTfacets, " - center: ", facet); +#if qh_MAXoutside + if (facet->maxoutside > qh DISTround) + fprintf(fp, " - maxoutside: %10.7g\n", facet->maxoutside); +#endif + if (!SETempty_(facet->outsideset)) { + furthest= (pointT*)qh_setlast(facet->outsideset); + if (qh_setsize (facet->outsideset) < 6) { + fprintf(fp, " - outside set (furthest p%d):\n", qh_pointid(furthest)); + FOREACHpoint_(facet->outsideset) + qh_printpoint(fp, " ", point); + }else if (qh_setsize (facet->outsideset) < 21) { + qh_printpoints(fp, " - outside set:", facet->outsideset); + }else { + fprintf(fp, " - outside set: %d points.", qh_setsize(facet->outsideset)); + qh_printpoint(fp, " Furthest", furthest); + } +#if !qh_COMPUTEfurthest + fprintf(fp, " - furthest distance= %2.2g\n", facet->furthestdist); +#endif + } + if (!SETempty_(facet->coplanarset)) { + furthest= (pointT*)qh_setlast(facet->coplanarset); + if (qh_setsize (facet->coplanarset) < 6) { + fprintf(fp, " - coplanar set (furthest p%d):\n", qh_pointid(furthest)); + FOREACHpoint_(facet->coplanarset) + qh_printpoint(fp, " ", point); + }else if (qh_setsize (facet->coplanarset) < 21) { + qh_printpoints(fp, " - coplanar set:", facet->coplanarset); + }else { + fprintf(fp, " - coplanar set: %d points.", qh_setsize(facet->coplanarset)); + qh_printpoint(fp, " Furthest", furthest); + } + zinc_(Zdistio); + qh_distplane (furthest, facet, &dist); + fprintf(fp, " furthest distance= %2.2g\n", dist); + } + qh_printvertices (fp, " - vertices:", facet->vertices); + fprintf(fp, " - neighboring facets: "); + FOREACHneighbor_(facet) { + if (neighbor == qh_MERGEridge) + fprintf(fp, " MERGE"); + else if (neighbor == qh_DUPLICATEridge) + fprintf(fp, " DUP"); + else + fprintf(fp, " f%d", neighbor->id); + } + fprintf(fp, "\n"); + qh RANDOMdist= qh old_randomdist; +} /* printfacetheader */ + + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printfacetridges">-</a> + + qh_printfacetridges( fp, facet ) + prints ridges of a facet to fp + + notes: + ridges printed in neighbor order + assumes the ridges exist + for 'f' output +*/ +void qh_printfacetridges(FILE *fp, facetT *facet) { + facetT *neighbor, **neighborp; + ridgeT *ridge, **ridgep; + int numridges= 0; + + + if (facet->visible && qh NEWfacets) { + fprintf(fp, " - ridges (ids may be garbage):"); + FOREACHridge_(facet->ridges) + fprintf(fp, " r%d", ridge->id); + fprintf(fp, "\n"); + }else { + fprintf(fp, " - ridges:\n"); + FOREACHridge_(facet->ridges) + ridge->seen= False; + if (qh hull_dim == 3) { + ridge= SETfirstt_(facet->ridges, ridgeT); + while (ridge && !ridge->seen) { + ridge->seen= True; + qh_printridge(fp, ridge); + numridges++; + ridge= qh_nextridge3d (ridge, facet, NULL); + } + }else { + FOREACHneighbor_(facet) { + FOREACHridge_(facet->ridges) { + if (otherfacet_(ridge,facet) == neighbor) { + ridge->seen= True; + qh_printridge(fp, ridge); + numridges++; + } + } + } + } + if (numridges != qh_setsize (facet->ridges)) { + fprintf (fp, " - all ridges:"); + FOREACHridge_(facet->ridges) + fprintf (fp, " r%d", ridge->id); + fprintf (fp, "\n"); + } + FOREACHridge_(facet->ridges) { + if (!ridge->seen) + qh_printridge(fp, ridge); + } + } +} /* printfacetridges */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printfacets">-</a> + + qh_printfacets( fp, format, facetlist, facets, printall ) + prints facetlist and/or facet set in output format + + notes: + also used for specialized formats ('FO' and summary) + turns off 'Rn' option since want actual numbers +*/ +void qh_printfacets(FILE *fp, int format, facetT *facetlist, setT *facets, boolT printall) { + int numfacets, numsimplicial, numridges, totneighbors, numcoplanars, numtricoplanars; + facetT *facet, **facetp; + setT *vertices; + coordT *center; + realT outerplane, innerplane; + + qh old_randomdist= qh RANDOMdist; + qh RANDOMdist= False; + if (qh CDDoutput && (format == qh_PRINTcentrums || format == qh_PRINTpointintersect || format == qh_PRINToff)) + fprintf (qh ferr, "qhull warning: CDD format is not available for centrums, halfspace\nintersections, and OFF file format.\n"); + if (format == qh_PRINTnone) + ; /* print nothing */ + else if (format == qh_PRINTaverage) { + vertices= qh_facetvertices (facetlist, facets, printall); + center= qh_getcenter (vertices); + fprintf (fp, "%d 1\n", qh hull_dim); + qh_printpointid (fp, NULL, qh hull_dim, center, -1); + qh_memfree (center, qh normal_size); + qh_settempfree (&vertices); + }else if (format == qh_PRINTextremes) { + if (qh DELAUNAY) + qh_printextremes_d (fp, facetlist, facets, printall); + else if (qh hull_dim == 2) + qh_printextremes_2d (fp, facetlist, facets, printall); + else + qh_printextremes (fp, facetlist, facets, printall); + }else if (format == qh_PRINToptions) + fprintf(fp, "Options selected for Qhull %s:\n%s\n", qh_VERSION, qh qhull_options); + else if (format == qh_PRINTpoints && !qh VORONOI) + qh_printpoints_out (fp, facetlist, facets, printall); + else if (format == qh_PRINTqhull) + fprintf (fp, "%s | %s\n", qh rbox_command, qh qhull_command); + else if (format == qh_PRINTsize) { + fprintf (fp, "0\n2 "); + fprintf (fp, qh_REAL_1, qh totarea); + fprintf (fp, qh_REAL_1, qh totvol); + fprintf (fp, "\n"); + }else if (format == qh_PRINTsummary) { + qh_countfacets (facetlist, facets, printall, &numfacets, &numsimplicial, + &totneighbors, &numridges, &numcoplanars, &numtricoplanars); + vertices= qh_facetvertices (facetlist, facets, printall); + fprintf (fp, "10 %d %d %d %d %d %d %d %d %d %d\n2 ", qh hull_dim, + qh num_points + qh_setsize (qh other_points), + qh num_vertices, qh num_facets - qh num_visible, + qh_setsize (vertices), numfacets, numcoplanars, + numfacets - numsimplicial, zzval_(Zdelvertextot), + numtricoplanars); + qh_settempfree (&vertices); + qh_outerinner (NULL, &outerplane, &innerplane); + fprintf (fp, qh_REAL_2n, outerplane, innerplane); + }else if (format == qh_PRINTvneighbors) + qh_printvneighbors (fp, facetlist, facets, printall); + else if (qh VORONOI && format == qh_PRINToff) + qh_printvoronoi (fp, format, facetlist, facets, printall); + else if (qh VORONOI && format == qh_PRINTgeom) { + qh_printbegin (fp, format, facetlist, facets, printall); + qh_printvoronoi (fp, format, facetlist, facets, printall); + qh_printend (fp, format, facetlist, facets, printall); + }else if (qh VORONOI + && (format == qh_PRINTvertices || format == qh_PRINTinner || format == qh_PRINTouter)) + qh_printvdiagram (fp, format, facetlist, facets, printall); + else { + qh_printbegin (fp, format, facetlist, facets, printall); + FORALLfacet_(facetlist) + qh_printafacet (fp, format, facet, printall); + FOREACHfacet_(facets) + qh_printafacet (fp, format, facet, printall); + qh_printend (fp, format, facetlist, facets, printall); + } + qh RANDOMdist= qh old_randomdist; +} /* printfacets */ + + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printhelp_degenerate">-</a> + + qh_printhelp_degenerate( fp ) + prints descriptive message for precision error + + notes: + no message if qh_QUICKhelp +*/ +void qh_printhelp_degenerate(FILE *fp) { + + if (qh MERGEexact || qh PREmerge || qh JOGGLEmax < REALmax/2) + fprintf(fp, "\n\ +A Qhull error has occurred. Qhull should have corrected the above\n\ +precision error. Please send the input and all of the output to\n\ +qhull_bug@geom.umn.edu\n"); + else if (!qh_QUICKhelp) { + fprintf(fp, "\n\ +Precision problems were detected during construction of the convex hull.\n\ +This occurs because convex hull algorithms assume that calculations are\n\ +exact, but floating-point arithmetic has roundoff errors.\n\ +\n\ +To correct for precision problems, do not use 'Q0'. By default, Qhull\n\ +selects 'C-0' or 'Qx' and merges non-convex facets. With option 'QJ',\n\ +Qhull joggles the input to prevent precision problems. See \"Imprecision\n\ +in Qhull\" (qh-impre.htm).\n\ +\n\ +If you use 'Q0', the output may include\n\ +coplanar ridges, concave ridges, and flipped facets. In 4-d and higher,\n\ +Qhull may produce a ridge with four neighbors or two facets with the same \n\ +vertices. Qhull reports these events when they occur. It stops when a\n\ +concave ridge, flipped facet, or duplicate facet occurs.\n"); +#if REALfloat + fprintf (fp, "\ +\n\ +Qhull is currently using single precision arithmetic. The following\n\ +will probably remove the precision problems:\n\ + - recompile qhull for double precision (#define REALfloat 0 in user.h).\n"); +#endif + if (qh DELAUNAY && !qh SCALElast && qh MAXabs_coord > 1e4) + fprintf( fp, "\ +\n\ +When computing the Delaunay triangulation of coordinates > 1.0,\n\ + - use 'Qbb' to scale the last coordinate to [0,m] (max previous coordinate)\n"); + if (qh DELAUNAY && !qh ATinfinity) + fprintf( fp, "\ +When computing the Delaunay triangulation:\n\ + - use 'Qz' to add a point at-infinity. This reduces precision problems.\n"); + + fprintf(fp, "\ +\n\ +If you need triangular output:\n\ + - use option 'Qt' to triangulate the output\n\ + - use option 'QJ' to joggle the input points and remove precision errors\n\ + - use option 'Ft'. It triangulates non-simplicial facets with added points.\n\ +\n\ +If you must use 'Q0',\n\ +try one or more of the following options. They can not guarantee an output.\n\ + - use 'QbB' to scale the input to a cube.\n\ + - use 'Po' to produce output and prevent partitioning for flipped facets\n\ + - use 'V0' to set min. distance to visible facet as 0 instead of roundoff\n\ + - use 'En' to specify a maximum roundoff error less than %2.2g.\n\ + - options 'Qf', 'Qbb', and 'QR0' may also help\n", + qh DISTround); + fprintf(fp, "\ +\n\ +To guarantee simplicial output:\n\ + - use option 'Qt' to triangulate the output\n\ + - use option 'QJ' to joggle the input points and remove precision errors\n\ + - use option 'Ft' to triangulate the output by adding points\n\ + - use exact arithmetic (see \"Imprecision in Qhull\", qh-impre.htm)\n\ +"); + } +} /* printhelp_degenerate */ + + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printhelp_singular">-</a> + + qh_printhelp_singular( fp ) + prints descriptive message for singular input +*/ +void qh_printhelp_singular(FILE *fp) { + facetT *facet; + vertexT *vertex, **vertexp; + realT min, max, *coord, dist; + int i,k; + + fprintf(fp, "\n\ +The input to qhull appears to be less than %d dimensional, or a\n\ +computation has overflowed.\n\n\ +Qhull could not construct a clearly convex simplex from points:\n", + qh hull_dim); + qh_printvertexlist (fp, "", qh facet_list, NULL, qh_ALL); + if (!qh_QUICKhelp) + fprintf(fp, "\n\ +The center point is coplanar with a facet, or a vertex is coplanar\n\ +with a neighboring facet. The maximum round off error for\n\ +computing distances is %2.2g. The center point, facets and distances\n\ +to the center point are as follows:\n\n", qh DISTround); + qh_printpointid (fp, "center point", qh hull_dim, qh interior_point, -1); + fprintf (fp, "\n"); + FORALLfacets { + fprintf (fp, "facet"); + FOREACHvertex_(facet->vertices) + fprintf (fp, " p%d", qh_pointid(vertex->point)); + zinc_(Zdistio); + qh_distplane(qh interior_point, facet, &dist); + fprintf (fp, " distance= %4.2g\n", dist); + } + if (!qh_QUICKhelp) { + if (qh HALFspace) + fprintf (fp, "\n\ +These points are the dual of the given halfspaces. They indicate that\n\ +the intersection is degenerate.\n"); + fprintf (fp,"\n\ +These points either have a maximum or minimum x-coordinate, or\n\ +they maximize the determinant for k coordinates. Trial points\n\ +are first selected from points that maximize a coordinate.\n"); + if (qh hull_dim >= qh_INITIALmax) + fprintf (fp, "\n\ +Because of the high dimension, the min x-coordinate and max-coordinate\n\ +points are used if the determinant is non-zero. Option 'Qs' will\n\ +do a better, though much slower, job. Instead of 'Qs', you can change\n\ +the points by randomly rotating the input with 'QR0'.\n"); + } + fprintf (fp, "\nThe min and max coordinates for each dimension are:\n"); + for (k=0; k < qh hull_dim; k++) { + min= REALmax; + max= -REALmin; + for (i=qh num_points, coord= qh first_point+k; i--; coord += qh hull_dim) { + maximize_(max, *coord); + minimize_(min, *coord); + } + fprintf (fp, " %d: %8.4g %8.4g difference= %4.4g\n", k, min, max, max-min); + } + if (!qh_QUICKhelp) { + fprintf (fp, "\n\ +If the input should be full dimensional, you have several options that\n\ +may determine an initial simplex:\n\ + - use 'QJ' to joggle the input and make it full dimensional\n\ + - use 'QbB' to scale the points to the unit cube\n\ + - use 'QR0' to randomly rotate the input for different maximum points\n\ + - use 'Qs' to search all points for the initial simplex\n\ + - use 'En' to specify a maximum roundoff error less than %2.2g.\n\ + - trace execution with 'T3' to see the determinant for each point.\n", + qh DISTround); +#if REALfloat + fprintf (fp, "\ + - recompile qhull for double precision (#define REALfloat 0 in qhull.h).\n"); +#endif + fprintf (fp, "\n\ +If the input is lower dimensional:\n\ + - use 'QJ' to joggle the input and make it full dimensional\n\ + - use 'Qbk:0Bk:0' to delete coordinate k from the input. You should\n\ + pick the coordinate with the least range. The hull will have the\n\ + correct topology.\n\ + - determine the flat containing the points, rotate the points\n\ + into a coordinate plane, and delete the other coordinates.\n\ + - add one or more points to make the input full dimensional.\n\ +"); + if (qh DELAUNAY && !qh ATinfinity) + fprintf (fp, "\n\n\ +This is a Delaunay triangulation and the input is co-circular or co-spherical:\n\ + - use 'Qz' to add a point \"at infinity\" (i.e., above the paraboloid)\n\ + - or use 'QJ' to joggle the input and avoid co-circular data\n"); + } +} /* printhelp_singular */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printhyperplaneintersection">-</a> + + qh_printhyperplaneintersection( fp, facet1, facet2, vertices, color ) + print Geomview OFF or 4OFF for the intersection of two hyperplanes in 3-d or 4-d +*/ +void qh_printhyperplaneintersection(FILE *fp, facetT *facet1, facetT *facet2, + setT *vertices, realT color[3]) { + realT costheta, denominator, dist1, dist2, s, t, mindenom, p[4]; + vertexT *vertex, **vertexp; + int i, k; + boolT nearzero1, nearzero2; + + costheta= qh_getangle(facet1->normal, facet2->normal); + denominator= 1 - costheta * costheta; + i= qh_setsize(vertices); + if (qh hull_dim == 3) + fprintf(fp, "VECT 1 %d 1 %d 1 ", i, i); + else if (qh hull_dim == 4 && qh DROPdim >= 0) + fprintf(fp, "OFF 3 1 1 "); + else + qh printoutvar++; + fprintf (fp, "# intersect f%d f%d\n", facet1->id, facet2->id); + mindenom= 1 / (10.0 * qh MAXabs_coord); + FOREACHvertex_(vertices) { + zadd_(Zdistio, 2); + qh_distplane(vertex->point, facet1, &dist1); + qh_distplane(vertex->point, facet2, &dist2); + s= qh_divzero (-dist1 + costheta * dist2, denominator,mindenom,&nearzero1); + t= qh_divzero (-dist2 + costheta * dist1, denominator,mindenom,&nearzero2); + if (nearzero1 || nearzero2) + s= t= 0.0; + for(k= qh hull_dim; k--; ) + p[k]= vertex->point[k] + facet1->normal[k] * s + facet2->normal[k] * t; + if (qh PRINTdim <= 3) { + qh_projectdim3 (p, p); + fprintf(fp, "%8.4g %8.4g %8.4g # ", p[0], p[1], p[2]); + }else + fprintf(fp, "%8.4g %8.4g %8.4g %8.4g # ", p[0], p[1], p[2], p[3]); + if (nearzero1+nearzero2) + fprintf (fp, "p%d (coplanar facets)\n", qh_pointid (vertex->point)); + else + fprintf (fp, "projected p%d\n", qh_pointid (vertex->point)); + } + if (qh hull_dim == 3) + fprintf(fp, "%8.4g %8.4g %8.4g 1.0\n", color[0], color[1], color[2]); + else if (qh hull_dim == 4 && qh DROPdim >= 0) + fprintf(fp, "3 0 1 2 %8.4g %8.4g %8.4g 1.0\n", color[0], color[1], color[2]); +} /* printhyperplaneintersection */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printline3geom">-</a> + + qh_printline3geom( fp, pointA, pointB, color ) + prints a line as a VECT + prints 0's for qh.DROPdim + + notes: + if pointA == pointB, + it's a 1 point VECT +*/ +void qh_printline3geom (FILE *fp, pointT *pointA, pointT *pointB, realT color[3]) { + int k; + realT pA[4], pB[4]; + + qh_projectdim3(pointA, pA); + qh_projectdim3(pointB, pB); + if ((fabs(pA[0] - pB[0]) > 1e-3) || + (fabs(pA[1] - pB[1]) > 1e-3) || + (fabs(pA[2] - pB[2]) > 1e-3)) { + fprintf (fp, "VECT 1 2 1 2 1\n"); + for (k= 0; k < 3; k++) + fprintf (fp, "%8.4g ", pB[k]); + fprintf (fp, " # p%d\n", qh_pointid (pointB)); + }else + fprintf (fp, "VECT 1 1 1 1 1\n"); + for (k=0; k < 3; k++) + fprintf (fp, "%8.4g ", pA[k]); + fprintf (fp, " # p%d\n", qh_pointid (pointA)); + fprintf (fp, "%8.4g %8.4g %8.4g 1\n", color[0], color[1], color[2]); +} + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printneighborhood">-</a> + + qh_printneighborhood( fp, format, facetA, facetB, printall ) + print neighborhood of one or two facets + + notes: + calls qh_findgood_all() + bumps qh.visit_id +*/ +void qh_printneighborhood (FILE *fp, int format, facetT *facetA, facetT *facetB, boolT printall) { + facetT *neighbor, **neighborp, *facet; + setT *facets; + + if (format == qh_PRINTnone) + return; + qh_findgood_all (qh facet_list); + if (facetA == facetB) + facetB= NULL; + facets= qh_settemp (2*(qh_setsize (facetA->neighbors)+1)); + qh visit_id++; + for (facet= facetA; facet; facet= ((facet == facetA) ? facetB : NULL)) { + if (facet->visitid != qh visit_id) { + facet->visitid= qh visit_id; + qh_setappend (&facets, facet); + } + FOREACHneighbor_(facet) { + if (neighbor->visitid == qh visit_id) + continue; + neighbor->visitid= qh visit_id; + if (printall || !qh_skipfacet (neighbor)) + qh_setappend (&facets, neighbor); + } + } + qh_printfacets (fp, format, NULL, facets, printall); + qh_settempfree (&facets); +} /* printneighborhood */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printpoint">-</a> + + qh_printpoint( fp, string, point ) + qh_printpointid( fp, string, dim, point, id ) + prints the coordinates of a point + + returns: + if string is defined + prints 'string p%d' (skips p%d if id=-1) + + notes: + nop if point is NULL + prints id unless it is undefined (-1) +*/ +void qh_printpoint(FILE *fp, char *string, pointT *point) { + int id= qh_pointid( point); + + qh_printpointid( fp, string, qh hull_dim, point, id); +} /* printpoint */ + +void qh_printpointid(FILE *fp, char *string, int dim, pointT *point, int id) { + int k; + realT r; /*bug fix*/ + + if (!point) + return; + if (string) { + fputs (string, fp); + if (id != -1) + fprintf(fp, " p%d: ", id); + } + for(k= dim; k--; ) { + r= *point++; + if (string) + fprintf(fp, " %8.4g", r); + else + fprintf(fp, qh_REAL_1, r); + } + fprintf(fp, "\n"); +} /* printpointid */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printpoint3">-</a> + + qh_printpoint3( fp, point ) + prints 2-d, 3-d, or 4-d point as Geomview 3-d coordinates +*/ +void qh_printpoint3 (FILE *fp, pointT *point) { + int k; + realT p[4]; + + qh_projectdim3 (point, p); + for (k=0; k < 3; k++) + fprintf (fp, "%8.4g ", p[k]); + fprintf (fp, " # p%d\n", qh_pointid (point)); +} /* printpoint3 */ + +/*---------------------------------------- +-printpoints- print pointids for a set of points starting at index + see geom.c +*/ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printpoints_out">-</a> + + qh_printpoints_out( fp, facetlist, facets, printall ) + prints vertices, coplanar/inside points, for facets by their point coordinates + allows qh.CDDoutput + + notes: + same format as qhull input + if no coplanar/interior points, + same order as qh_printextremes +*/ +void qh_printpoints_out (FILE *fp, facetT *facetlist, setT *facets, int printall) { + int allpoints= qh num_points + qh_setsize (qh other_points); + int numpoints=0, point_i, point_n; + setT *vertices, *points; + facetT *facet, **facetp; + pointT *point, **pointp; + vertexT *vertex, **vertexp; + int id; + + points= qh_settemp (allpoints); + qh_setzero (points, 0, allpoints); + vertices= qh_facetvertices (facetlist, facets, printall); + FOREACHvertex_(vertices) { + id= qh_pointid (vertex->point); + if (id >= 0) + SETelem_(points, id)= vertex->point; + } + if (qh KEEPinside || qh KEEPcoplanar || qh KEEPnearinside) { + FORALLfacet_(facetlist) { + if (!printall && qh_skipfacet(facet)) + continue; + FOREACHpoint_(facet->coplanarset) { + id= qh_pointid (point); + if (id >= 0) + SETelem_(points, id)= point; + } + } + FOREACHfacet_(facets) { + if (!printall && qh_skipfacet(facet)) + continue; + FOREACHpoint_(facet->coplanarset) { + id= qh_pointid (point); + if (id >= 0) + SETelem_(points, id)= point; + } + } + } + qh_settempfree (&vertices); + FOREACHpoint_i_(points) { + if (point) + numpoints++; + } + if (qh CDDoutput) + fprintf (fp, "%s | %s\nbegin\n%d %d real\n", qh rbox_command, + qh qhull_command, numpoints, qh hull_dim + 1); + else + fprintf (fp, "%d\n%d\n", qh hull_dim, numpoints); + FOREACHpoint_i_(points) { + if (point) { + if (qh CDDoutput) + fprintf (fp, "1 "); + qh_printpoint (fp, NULL, point); + } + } + if (qh CDDoutput) + fprintf (fp, "end\n"); + qh_settempfree (&points); +} /* printpoints_out */ + + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printpointvect">-</a> + + qh_printpointvect( fp, point, normal, center, radius, color ) + prints a 2-d, 3-d, or 4-d point as 3-d VECT's relative to normal or to center point +*/ +void qh_printpointvect (FILE *fp, pointT *point, coordT *normal, pointT *center, realT radius, realT color[3]) { + realT diff[4], pointA[4]; + int k; + + for (k= qh hull_dim; k--; ) { + if (center) + diff[k]= point[k]-center[k]; + else if (normal) + diff[k]= normal[k]; + else + diff[k]= 0; + } + if (center) + qh_normalize2 (diff, qh hull_dim, True, NULL, NULL); + for (k= qh hull_dim; k--; ) + pointA[k]= point[k]+diff[k] * radius; + qh_printline3geom (fp, point, pointA, color); +} /* printpointvect */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printpointvect2">-</a> + + qh_printpointvect2( fp, point, normal, center, radius ) + prints a 2-d, 3-d, or 4-d point as 2 3-d VECT's for an imprecise point +*/ +void qh_printpointvect2 (FILE *fp, pointT *point, coordT *normal, pointT *center, realT radius) { + realT red[3]={1, 0, 0}, yellow[3]={1, 1, 0}; + + qh_printpointvect (fp, point, normal, center, radius, red); + qh_printpointvect (fp, point, normal, center, -radius, yellow); +} /* printpointvect2 */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printridge">-</a> + + qh_printridge( fp, ridge ) + prints the information in a ridge + + notes: + for qh_printfacetridges() +*/ +void qh_printridge(FILE *fp, ridgeT *ridge) { + + fprintf(fp, " - r%d", ridge->id); + if (ridge->tested) + fprintf (fp, " tested"); + if (ridge->nonconvex) + fprintf (fp, " nonconvex"); + fprintf (fp, "\n"); + qh_printvertices (fp, " vertices:", ridge->vertices); + if (ridge->top && ridge->bottom) + fprintf(fp, " between f%d and f%d\n", + ridge->top->id, ridge->bottom->id); +} /* printridge */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printspheres">-</a> + + qh_printspheres( fp, vertices, radius ) + prints 3-d vertices as OFF spheres + + notes: + inflated octahedron from Stuart Levy earth/mksphere2 +*/ +void qh_printspheres(FILE *fp, setT *vertices, realT radius) { + vertexT *vertex, **vertexp; + + qh printoutnum++; + fprintf (fp, "{appearance {-edge -normal normscale 0} {\n\ +INST geom {define vsphere OFF\n\ +18 32 48\n\ +\n\ +0 0 1\n\ +1 0 0\n\ +0 1 0\n\ +-1 0 0\n\ +0 -1 0\n\ +0 0 -1\n\ +0.707107 0 0.707107\n\ +0 -0.707107 0.707107\n\ +0.707107 -0.707107 0\n\ +-0.707107 0 0.707107\n\ +-0.707107 -0.707107 0\n\ +0 0.707107 0.707107\n\ +-0.707107 0.707107 0\n\ +0.707107 0.707107 0\n\ +0.707107 0 -0.707107\n\ +0 0.707107 -0.707107\n\ +-0.707107 0 -0.707107\n\ +0 -0.707107 -0.707107\n\ +\n\ +3 0 6 11\n\ +3 0 7 6 \n\ +3 0 9 7 \n\ +3 0 11 9\n\ +3 1 6 8 \n\ +3 1 8 14\n\ +3 1 13 6\n\ +3 1 14 13\n\ +3 2 11 13\n\ +3 2 12 11\n\ +3 2 13 15\n\ +3 2 15 12\n\ +3 3 9 12\n\ +3 3 10 9\n\ +3 3 12 16\n\ +3 3 16 10\n\ +3 4 7 10\n\ +3 4 8 7\n\ +3 4 10 17\n\ +3 4 17 8\n\ +3 5 14 17\n\ +3 5 15 14\n\ +3 5 16 15\n\ +3 5 17 16\n\ +3 6 13 11\n\ +3 7 8 6\n\ +3 9 10 7\n\ +3 11 12 9\n\ +3 14 8 17\n\ +3 15 13 14\n\ +3 16 12 15\n\ +3 17 10 16\n} transforms { TLIST\n"); + FOREACHvertex_(vertices) { + fprintf(fp, "%8.4g 0 0 0 # v%d\n 0 %8.4g 0 0\n0 0 %8.4g 0\n", + radius, vertex->id, radius, radius); + qh_printpoint3 (fp, vertex->point); + fprintf (fp, "1\n"); + } + fprintf (fp, "}}}\n"); +} /* printspheres */ + + +/*---------------------------------------------- +-printsummary- + see qhull.c +*/ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printvdiagram">-</a> + + qh_printvdiagram( fp, format, facetlist, facets, printall ) + print voronoi diagram + # of pairs of input sites + #indices site1 site2 vertex1 ... + + sites indexed by input point id + point 0 is the first input point + vertices indexed by 'o' and 'p' order + vertex 0 is the 'vertex-at-infinity' + vertex 1 is the first Voronoi vertex + + see: + qh_printvoronoi() + qh_eachvoronoi_all() + + notes: + if all facets are upperdelaunay, + prints upper hull (furthest-site Voronoi diagram) +*/ +void qh_printvdiagram (FILE *fp, int format, facetT *facetlist, setT *facets, boolT printall) { + setT *vertices; + int totcount, numcenters; + boolT islower; + qh_RIDGE innerouter= qh_RIDGEall; + printvridgeT printvridge= NULL; + + if (format == qh_PRINTvertices) { + innerouter= qh_RIDGEall; + printvridge= qh_printvridge; + }else if (format == qh_PRINTinner) { + innerouter= qh_RIDGEinner; + printvridge= qh_printvnorm; + }else if (format == qh_PRINTouter) { + innerouter= qh_RIDGEouter; + printvridge= qh_printvnorm; + }else { + fprintf(qh ferr, "qh_printvdiagram: unknown print format %d.\n", format); + qh_errexit (qh_ERRinput, NULL, NULL); + } + vertices= qh_markvoronoi (facetlist, facets, printall, &islower, &numcenters); + totcount= qh_printvdiagram2 (NULL, NULL, vertices, innerouter, False); + fprintf (fp, "%d\n", totcount); + totcount= qh_printvdiagram2 (fp, printvridge, vertices, innerouter, True /* inorder*/); + qh_settempfree (&vertices); +#if 0 /* for testing qh_eachvoronoi_all */ + fprintf (fp, "\n"); + totcount= qh_eachvoronoi_all(fp, printvridge, qh UPPERdelaunay, innerouter, True /* inorder*/); + fprintf (fp, "%d\n", totcount); +#endif +} /* printvdiagram */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printvdiagram2">-</a> + + qh_printvdiagram2( fp, printvridge, vertices, innerouter, inorder ) + visit all pairs of input sites (vertices) for selected Voronoi vertices + vertices may include NULLs + + innerouter: + qh_RIDGEall print inner ridges (bounded) and outer ridges (unbounded) + qh_RIDGEinner print only inner ridges + qh_RIDGEouter print only outer ridges + + inorder: + print 3-d Voronoi vertices in order + + assumes: + qh_markvoronoi marked facet->visitid for Voronoi vertices + all facet->seen= False + all facet->seen2= True + + returns: + total number of Voronoi ridges + if printvridge, + calls printvridge( fp, vertex, vertexA, centers) for each ridge + [see qh_eachvoronoi()] + + see: + qh_eachvoronoi_all() +*/ +int qh_printvdiagram2 (FILE *fp, printvridgeT printvridge, setT *vertices, qh_RIDGE innerouter, boolT inorder) { + int totcount= 0; + int vertex_i, vertex_n; + vertexT *vertex; + + FORALLvertices + vertex->seen= False; + FOREACHvertex_i_(vertices) { + if (vertex) { + if (qh GOODvertex > 0 && qh_pointid(vertex->point)+1 != qh GOODvertex) + continue; + totcount += qh_eachvoronoi (fp, printvridge, vertex, !qh_ALL, innerouter, inorder); + } + } + return totcount; +} /* printvdiagram2 */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printvertex">-</a> + + qh_printvertex( fp, vertex ) + prints the information in a vertex +*/ +void qh_printvertex(FILE *fp, vertexT *vertex) { + pointT *point; + int k, count= 0; + facetT *neighbor, **neighborp; + realT r; /*bug fix*/ + + if (!vertex) { + fprintf (fp, " NULLvertex\n"); + return; + } + fprintf(fp, "- p%d (v%d):", qh_pointid(vertex->point), vertex->id); + point= vertex->point; + if (point) { + for(k= qh hull_dim; k--; ) { + r= *point++; + fprintf(fp, " %5.2g", r); + } + } + if (vertex->deleted) + fprintf(fp, " deleted"); + if (vertex->delridge) + fprintf (fp, " ridgedeleted"); + fprintf(fp, "\n"); + if (vertex->neighbors) { + fprintf(fp, " neighbors:"); + FOREACHneighbor_(vertex) { + if (++count % 100 == 0) + fprintf (fp, "\n "); + fprintf(fp, " f%d", neighbor->id); + } + fprintf(fp, "\n"); + } +} /* printvertex */ + + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printvertexlist">-</a> + + qh_printvertexlist( fp, string, facetlist, facets, printall ) + prints vertices used by a facetlist or facet set + tests qh_skipfacet() if !printall +*/ +void qh_printvertexlist (FILE *fp, char* string, facetT *facetlist, + setT *facets, boolT printall) { + vertexT *vertex, **vertexp; + setT *vertices; + + vertices= qh_facetvertices (facetlist, facets, printall); + fputs (string, fp); + FOREACHvertex_(vertices) + qh_printvertex(fp, vertex); + qh_settempfree (&vertices); +} /* printvertexlist */ + + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printvertices">-</a> + + qh_printvertices( fp, string, vertices ) + prints vertices in a set +*/ +void qh_printvertices(FILE *fp, char* string, setT *vertices) { + vertexT *vertex, **vertexp; + + fputs (string, fp); + FOREACHvertex_(vertices) + fprintf (fp, " p%d (v%d)", qh_pointid(vertex->point), vertex->id); + fprintf(fp, "\n"); +} /* printvertices */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printvneighbors">-</a> + + qh_printvneighbors( fp, facetlist, facets, printall ) + print vertex neighbors of vertices in facetlist and facets ('FN') + + notes: + qh_countfacets clears facet->visitid for non-printed facets + + design: + collect facet count and related statistics + if necessary, build neighbor sets for each vertex + collect vertices in facetlist and facets + build a point array for point->vertex and point->coplanar facet + for each point + list vertex neighbors or coplanar facet +*/ +void qh_printvneighbors (FILE *fp, facetT* facetlist, setT *facets, boolT printall) { + int numfacets, numsimplicial, numridges, totneighbors, numneighbors, numcoplanars, numtricoplanars; + setT *vertices, *vertex_points, *coplanar_points; + int numpoints= qh num_points + qh_setsize (qh other_points); + vertexT *vertex, **vertexp; + int vertex_i, vertex_n; + facetT *facet, **facetp, *neighbor, **neighborp; + pointT *point, **pointp; + + qh_countfacets (facetlist, facets, printall, &numfacets, &numsimplicial, + &totneighbors, &numridges, &numcoplanars, &numtricoplanars); /* sets facet->visitid */ + fprintf (fp, "%d\n", numpoints); + qh_vertexneighbors(); + vertices= qh_facetvertices (facetlist, facets, printall); + vertex_points= qh_settemp (numpoints); + coplanar_points= qh_settemp (numpoints); + qh_setzero (vertex_points, 0, numpoints); + qh_setzero (coplanar_points, 0, numpoints); + FOREACHvertex_(vertices) + qh_point_add (vertex_points, vertex->point, vertex); + FORALLfacet_(facetlist) { + FOREACHpoint_(facet->coplanarset) + qh_point_add (coplanar_points, point, facet); + } + FOREACHfacet_(facets) { + FOREACHpoint_(facet->coplanarset) + qh_point_add (coplanar_points, point, facet); + } + FOREACHvertex_i_(vertex_points) { + if (vertex) { + numneighbors= qh_setsize (vertex->neighbors); + fprintf (fp, "%d", numneighbors); + if (qh hull_dim == 3) + qh_order_vertexneighbors (vertex); + else if (qh hull_dim >= 4) + qsort (SETaddr_(vertex->neighbors, facetT), numneighbors, + sizeof (facetT *), qh_compare_facetvisit); + FOREACHneighbor_(vertex) + fprintf (fp, " %d", + neighbor->visitid ? neighbor->visitid - 1 : - neighbor->id); + fprintf (fp, "\n"); + }else if ((facet= SETelemt_(coplanar_points, vertex_i, facetT))) + fprintf (fp, "1 %d\n", + facet->visitid ? facet->visitid - 1 : - facet->id); + else + fprintf (fp, "0\n"); + } + qh_settempfree (&coplanar_points); + qh_settempfree (&vertex_points); + qh_settempfree (&vertices); +} /* printvneighbors */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printvoronoi">-</a> + + qh_printvoronoi( fp, format, facetlist, facets, printall ) + print voronoi diagram in 'o' or 'G' format + for 'o' format + prints voronoi centers for each facet and for infinity + for each vertex, lists ids of printed facets or infinity + assumes facetlist and facets are disjoint + for 'G' format + prints an OFF object + adds a 0 coordinate to center + prints infinity but does not list in vertices + + see: + qh_printvdiagram() + + notes: + if 'o', + prints a line for each point except "at-infinity" + if all facets are upperdelaunay, + reverses lower and upper hull +*/ +void qh_printvoronoi (FILE *fp, int format, facetT *facetlist, setT *facets, boolT printall) { + int k, numcenters, numvertices= 0, numneighbors, numinf, vid=1, vertex_i, vertex_n; + facetT *facet, **facetp, *neighbor, **neighborp; + setT *vertices; + vertexT *vertex; + boolT islower; + unsigned int numfacets= (unsigned int) qh num_facets; + + vertices= qh_markvoronoi (facetlist, facets, printall, &islower, &numcenters); + FOREACHvertex_i_(vertices) { + if (vertex) { + numvertices++; + numneighbors = numinf = 0; + FOREACHneighbor_(vertex) { + if (neighbor->visitid == 0) + numinf= 1; + else if (neighbor->visitid < numfacets) + numneighbors++; + } + if (numinf && !numneighbors) { + SETelem_(vertices, vertex_i)= NULL; + numvertices--; + } + } + } + if (format == qh_PRINTgeom) + fprintf (fp, "{appearance {+edge -face} OFF %d %d 1 # Voronoi centers and cells\n", + numcenters, numvertices); + else + fprintf (fp, "%d\n%d %d 1\n", qh hull_dim-1, numcenters, qh_setsize(vertices)); + if (format == qh_PRINTgeom) { + for (k= qh hull_dim-1; k--; ) + fprintf (fp, qh_REAL_1, 0.0); + fprintf (fp, " 0 # infinity not used\n"); + }else { + for (k= qh hull_dim-1; k--; ) + fprintf (fp, qh_REAL_1, qh_INFINITE); + fprintf (fp, "\n"); + } + FORALLfacet_(facetlist) { + if (facet->visitid && facet->visitid < numfacets) { + if (format == qh_PRINTgeom) + fprintf (fp, "# %d f%d\n", vid++, facet->id); + qh_printcenter (fp, format, NULL, facet); + } + } + FOREACHfacet_(facets) { + if (facet->visitid && facet->visitid < numfacets) { + if (format == qh_PRINTgeom) + fprintf (fp, "# %d f%d\n", vid++, facet->id); + qh_printcenter (fp, format, NULL, facet); + } + } + FOREACHvertex_i_(vertices) { + numneighbors= 0; + numinf=0; + if (vertex) { + if (qh hull_dim == 3) + qh_order_vertexneighbors(vertex); + else if (qh hull_dim >= 4) + qsort (SETaddr_(vertex->neighbors, vertexT), + qh_setsize (vertex->neighbors), + sizeof (facetT *), qh_compare_facetvisit); + FOREACHneighbor_(vertex) { + if (neighbor->visitid == 0) + numinf= 1; + else if (neighbor->visitid < numfacets) + numneighbors++; + } + } + if (format == qh_PRINTgeom) { + if (vertex) { + fprintf (fp, "%d", numneighbors); + if (vertex) { + FOREACHneighbor_(vertex) { + if (neighbor->visitid && neighbor->visitid < numfacets) + fprintf (fp, " %d", neighbor->visitid); + } + } + fprintf (fp, " # p%d (v%d)\n", vertex_i, vertex->id); + }else + fprintf (fp, " # p%d is coplanar or isolated\n", vertex_i); + }else { + if (numinf) + numneighbors++; + fprintf (fp, "%d", numneighbors); + if (vertex) { + FOREACHneighbor_(vertex) { + if (neighbor->visitid == 0) { + if (numinf) { + numinf= 0; + fprintf (fp, " %d", neighbor->visitid); + } + }else if (neighbor->visitid < numfacets) + fprintf (fp, " %d", neighbor->visitid); + } + } + fprintf (fp, "\n"); + } + } + if (format == qh_PRINTgeom) + fprintf (fp, "}\n"); + qh_settempfree (&vertices); +} /* printvoronoi */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printvnorm">-</a> + + qh_printvnorm( fp, vertex, vertexA, centers, unbounded ) + print one separating plane of the Voronoi diagram for a pair of input sites + unbounded==True if centers includes vertex-at-infinity + + assumes: + qh_ASvoronoi and qh_vertexneighbors() already set + + see: + qh_printvdiagram() + qh_eachvoronoi() +*/ +void qh_printvnorm (FILE *fp, vertexT *vertex, vertexT *vertexA, setT *centers, boolT unbounded) { + pointT *normal; + realT offset; + int k; + + normal= qh_detvnorm (vertex, vertexA, centers, &offset); + fprintf (fp, "%d %d %d ", + 2+qh hull_dim, qh_pointid (vertex->point), qh_pointid (vertexA->point)); + for (k= 0; k< qh hull_dim-1; k++) + fprintf (fp, qh_REAL_1, normal[k]); + fprintf (fp, qh_REAL_1, offset); + fprintf (fp, "\n"); +} /* printvnorm */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="printvridge">-</a> + + qh_printvridge( fp, vertex, vertexA, centers, unbounded ) + print one ridge of the Voronoi diagram for a pair of input sites + unbounded==True if centers includes vertex-at-infinity + + see: + qh_printvdiagram() + + notes: + the user may use a different function +*/ +void qh_printvridge (FILE *fp, vertexT *vertex, vertexT *vertexA, setT *centers, boolT unbounded) { + facetT *facet, **facetp; + + fprintf (fp, "%d %d %d", qh_setsize (centers)+2, + qh_pointid (vertex->point), qh_pointid (vertexA->point)); + FOREACHfacet_(centers) + fprintf (fp, " %d", facet->visitid); + fprintf (fp, "\n"); +} /* printvridge */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="projectdim3">-</a> + + qh_projectdim3( source, destination ) + project 2-d 3-d or 4-d point to a 3-d point + uses qh.DROPdim and qh.hull_dim + source and destination may be the same + + notes: + allocate 4 elements to destination just in case +*/ +void qh_projectdim3 (pointT *source, pointT *destination) { + int i,k; + + for (k= 0, i=0; k < qh hull_dim; k++) { + if (qh hull_dim == 4) { + if (k != qh DROPdim) + destination[i++]= source[k]; + }else if (k == qh DROPdim) + destination[i++]= 0; + else + destination[i++]= source[k]; + } + while (i < 3) + destination[i++]= 0.0; +} /* projectdim3 */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="readfeasible">-</a> + + qh_readfeasible( dim, remainder ) + read feasible point from remainder string and qh.fin + + returns: + number of lines read from qh.fin + sets qh.FEASIBLEpoint with malloc'd coordinates + + notes: + checks for qh.HALFspace + assumes dim > 1 + + see: + qh_setfeasible +*/ +int qh_readfeasible (int dim, char *remainder) { + boolT isfirst= True; + int linecount= 0, tokcount= 0; + char *s, *t, firstline[qh_MAXfirst+1]; + coordT *coords, value; + + if (!qh HALFspace) { + fprintf (qh ferr, "qhull input error: feasible point (dim 1 coords) is only valid for halfspace intersection\n"); + qh_errexit (qh_ERRinput, NULL, NULL); + } + if (qh feasible_string) + fprintf (qh ferr, "qhull input warning: feasible point (dim 1 coords) overrides 'Hn,n,n' feasible point for halfspace intersection\n"); + if (!(qh feasible_point= (coordT*)malloc (dim* sizeof(coordT)))) { + fprintf(qh ferr, "qhull error: insufficient memory for feasible point\n"); + qh_errexit(qh_ERRmem, NULL, NULL); + } + coords= qh feasible_point; + while ((s= (isfirst ? remainder : fgets(firstline, qh_MAXfirst, qh fin)))) { + if (isfirst) + isfirst= False; + else + linecount++; + while (*s) { + while (isspace(*s)) + s++; + value= qh_strtod (s, &t); + if (s == t) + break; + s= t; + *(coords++)= value; + if (++tokcount == dim) { + while (isspace (*s)) + s++; + qh_strtod (s, &t); + if (s != t) { + fprintf (qh ferr, "qhull input error: coordinates for feasible point do not finish out the line: %s\n", + s); + qh_errexit (qh_ERRinput, NULL, NULL); + } + return linecount; + } + } + } + fprintf (qh ferr, "qhull input error: only %d coordinates. Could not read %d-d feasible point.\n", + tokcount, dim); + qh_errexit (qh_ERRinput, NULL, NULL); + return 0; +} /* readfeasible */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="readpoints">-</a> + + qh_readpoints( numpoints, dimension, ismalloc ) + read points from qh.fin into qh.first_point, qh.num_points + qh.fin is lines of coordinates, one per vertex, first line number of points + if 'rbox D4', + gives message + if qh.ATinfinity, + adds point-at-infinity for Delaunay triangulations + + returns: + number of points, array of point coordinates, dimension, ismalloc True + if qh.DELAUNAY & !qh.PROJECTinput, projects points to paraboloid + and clears qh.PROJECTdelaunay + if qh.HALFspace, reads optional feasible point, reads halfspaces, + converts to dual. + + for feasible point in "cdd format" in 3-d: + 3 1 + coordinates + comments + begin + n 4 real/integer + ... + end + + notes: + dimension will change in qh_initqhull_globals if qh.PROJECTinput + uses malloc() since qh_mem not initialized + FIXUP: this routine needs rewriting +*/ +coordT *qh_readpoints(int *numpoints, int *dimension, boolT *ismalloc) { + coordT *points, *coords, *infinity= NULL; + realT paraboloid, maxboloid= -REALmax, value; + realT *coordp= NULL, *offsetp= NULL, *normalp= NULL; + char *s, *t, firstline[qh_MAXfirst+1]; + int diminput=0, numinput=0, dimfeasible= 0, newnum, k, tempi; + int firsttext=0, firstshort=0, firstlong=0, firstpoint=0; + int tokcount= 0, linecount=0, maxcount, coordcount=0; + boolT islong, isfirst= True, wasbegin= False; + boolT isdelaunay= qh DELAUNAY && !qh PROJECTinput; + + if (qh CDDinput) { + while ((s= fgets(firstline, qh_MAXfirst, qh fin))) { + linecount++; + if (qh HALFspace && linecount == 1 && isdigit(*s)) { + dimfeasible= qh_strtol (s, &s); + while (isspace(*s)) + s++; + if (qh_strtol (s, &s) == 1) + linecount += qh_readfeasible (dimfeasible, s); + else + dimfeasible= 0; + }else if (!memcmp (firstline, "begin", 5) || !memcmp (firstline, "BEGIN", 5)) + break; + else if (!*qh rbox_command) + strncat(qh rbox_command, s, sizeof (qh rbox_command)-1); + } + if (!s) { + fprintf (qh ferr, "qhull input error: missing \"begin\" for cdd-formated input\n"); + qh_errexit (qh_ERRinput, NULL, NULL); + } + } + while(!numinput && (s= fgets(firstline, qh_MAXfirst, qh fin))) { + linecount++; + if (!memcmp (s, "begin", 5) || !memcmp (s, "BEGIN", 5)) + wasbegin= True; + while (*s) { + while (isspace(*s)) + s++; + if (!*s) + break; + if (!isdigit(*s)) { + if (!*qh rbox_command) { + strncat(qh rbox_command, s, sizeof (qh rbox_command)-1); + firsttext= linecount; + } + break; + } + if (!diminput) + diminput= qh_strtol (s, &s); + else { + numinput= qh_strtol (s, &s); + if (numinput == 1 && diminput >= 2 && qh HALFspace && !qh CDDinput) { + linecount += qh_readfeasible (diminput, s); /* checks if ok */ + dimfeasible= diminput; + diminput= numinput= 0; + }else + break; + } + } + } + if (!s) { + fprintf(qh ferr, "qhull input error: short input file. Did not find dimension and number of points\n"); + qh_errexit(qh_ERRinput, NULL, NULL); + } + if (diminput > numinput) { + tempi= diminput; /* exchange dim and n, e.g., for cdd input format */ + diminput= numinput; + numinput= tempi; + } + if (diminput < 2) { + fprintf(qh ferr,"qhull input error: dimension %d (first number) should be at least 2\n", + diminput); + qh_errexit(qh_ERRinput, NULL, NULL); + } + if (isdelaunay) { + qh PROJECTdelaunay= False; + if (qh CDDinput) + *dimension= diminput; + else + *dimension= diminput+1; + *numpoints= numinput; + if (qh ATinfinity) + (*numpoints)++; + }else if (qh HALFspace) { + *dimension= diminput - 1; + *numpoints= numinput; + if (diminput < 3) { + fprintf(qh ferr,"qhull input error: dimension %d (first number, includes offset) should be at least 3 for halfspaces\n", + diminput); + qh_errexit(qh_ERRinput, NULL, NULL); + } + if (dimfeasible) { + if (dimfeasible != *dimension) { + fprintf(qh ferr,"qhull input error: dimension %d of feasible point is not one less than dimension %d for halfspaces\n", + dimfeasible, diminput); + qh_errexit(qh_ERRinput, NULL, NULL); + } + }else + qh_setfeasible (*dimension); + }else { + if (qh CDDinput) + *dimension= diminput-1; + else + *dimension= diminput; + *numpoints= numinput; + } + qh normal_size= *dimension * sizeof(coordT); /* for tracing with qh_printpoint */ + if (qh HALFspace) { + qh half_space= coordp= (coordT*) malloc (qh normal_size + sizeof(coordT)); + if (qh CDDinput) { + offsetp= qh half_space; + normalp= offsetp + 1; + }else { + normalp= qh half_space; + offsetp= normalp + *dimension; + } + } + qh maxline= diminput * (qh_REALdigits + 5); + maximize_(qh maxline, 500); + qh line= (char*)malloc ((qh maxline+1) * sizeof (char)); + *ismalloc= True; /* use malloc since memory not setup */ + coords= points= qh temp_malloc= + (coordT*)malloc((*numpoints)*(*dimension)*sizeof(coordT)); + if (!coords || !qh line || (qh HALFspace && !qh half_space)) { + fprintf(qh ferr, "qhull error: insufficient memory to read %d points\n", + numinput); + qh_errexit(qh_ERRmem, NULL, NULL); + } + if (isdelaunay && qh ATinfinity) { + infinity= points + numinput * (*dimension); + for (k= (*dimension) - 1; k--; ) + infinity[k]= 0.0; + } + maxcount= numinput * diminput; + paraboloid= 0.0; + while ((s= (isfirst ? s : fgets(qh line, qh maxline, qh fin)))) { + if (!isfirst) { + linecount++; + if (*s == 'e' || *s == 'E') { + if (!memcmp (s, "end", 3) || !memcmp (s, "END", 3)) { + if (qh CDDinput ) + break; + else if (wasbegin) + fprintf (qh ferr, "qhull input warning: the input appears to be in cdd format. If so, use 'Fd'\n"); + } + } + } + islong= False; + while (*s) { + while (isspace(*s)) + s++; + value= qh_strtod (s, &t); + if (s == t) { + if (!*qh rbox_command) + strncat(qh rbox_command, s, sizeof (qh rbox_command)-1); + if (*s && !firsttext) + firsttext= linecount; + if (!islong && !firstshort && coordcount) + firstshort= linecount; + break; + } + if (!firstpoint) + firstpoint= linecount; + s= t; + if (++tokcount > maxcount) + continue; + if (qh HALFspace) { + if (qh CDDinput) + *(coordp++)= -value; /* both coefficients and offset */ + else + *(coordp++)= value; + }else { + *(coords++)= value; + if (qh CDDinput && !coordcount) { + if (value != 1.0) { + fprintf (qh ferr, "qhull input error: for cdd format, point at line %d does not start with '1'\n", + linecount); + qh_errexit (qh_ERRinput, NULL, NULL); + } + coords--; + }else if (isdelaunay) { + paraboloid += value * value; + if (qh ATinfinity) { + if (qh CDDinput) + infinity[coordcount-1] += value; + else + infinity[coordcount] += value; + } + } + } + if (++coordcount == diminput) { + coordcount= 0; + if (isdelaunay) { + *(coords++)= paraboloid; + maximize_(maxboloid, paraboloid); + paraboloid= 0.0; + }else if (qh HALFspace) { + if (!qh_sethalfspace (*dimension, coords, &coords, normalp, offsetp, qh feasible_point)) { + fprintf (qh ferr, "The halfspace was on line %d\n", linecount); + if (wasbegin) + fprintf (qh ferr, "The input appears to be in cdd format. If so, you should use option 'Fd'\n"); + qh_errexit (qh_ERRinput, NULL, NULL); + } + coordp= qh half_space; + } + while (isspace(*s)) + s++; + if (*s) { + islong= True; + if (!firstlong) + firstlong= linecount; + } + } + } + if (!islong && !firstshort && coordcount) + firstshort= linecount; + if (!isfirst && s - qh line >= qh maxline) { + fprintf(qh ferr, "qhull input error: line %d contained more than %d characters\n", + linecount, (int) (s - qh line)); + qh_errexit(qh_ERRinput, NULL, NULL); + } + isfirst= False; + } + if (tokcount != maxcount) { + newnum= fmin_(numinput, tokcount/diminput); + fprintf(qh ferr,"\ +qhull warning: instead of %d %d-dimensional points, input contains\n\ +%d points and %d extra coordinates. Line %d is the first\npoint", + numinput, diminput, tokcount/diminput, tokcount % diminput, firstpoint); + if (firsttext) + fprintf(qh ferr, ", line %d is the first comment", firsttext); + if (firstshort) + fprintf(qh ferr, ", line %d is the first short\nline", firstshort); + if (firstlong) + fprintf(qh ferr, ", line %d is the first long line", firstlong); + fprintf(qh ferr, ". Continue with %d points.\n", newnum); + numinput= newnum; + if (isdelaunay && qh ATinfinity) { + for (k= tokcount % diminput; k--; ) + infinity[k] -= *(--coords); + *numpoints= newnum+1; + }else { + coords -= tokcount % diminput; + *numpoints= newnum; + } + } + if (isdelaunay && qh ATinfinity) { + for (k= (*dimension) -1; k--; ) + infinity[k] /= numinput; + if (coords == infinity) + coords += (*dimension) -1; + else { + for (k= 0; k < (*dimension) -1; k++) + *(coords++)= infinity[k]; + } + *(coords++)= maxboloid * 1.1; + } + if (qh rbox_command[0]) { + qh rbox_command[strlen(qh rbox_command)-1]= '\0'; + if (!strcmp (qh rbox_command, "./rbox D4")) + fprintf (qh ferr, "\n\ +This is the qhull test case. If any errors or core dumps occur,\n\ +recompile qhull with 'make new'. If errors still occur, there is\n\ +an incompatibility. You should try a different compiler. You can also\n\ +change the choices in user.h. If you discover the source of the problem,\n\ +please send mail to qhull_bug@geom.umn.edu.\n\ +\n\ +Type 'qhull' for a short list of options.\n"); + } + free (qh line); + qh line= NULL; + if (qh half_space) { + free (qh half_space); + qh half_space= NULL; + } + qh temp_malloc= NULL; + trace1((qh ferr,"qh_readpoints: read in %d %d-dimensional points\n", + numinput, diminput)); + return(points); +} /* readpoints */ + + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="setfeasible">-</a> + + qh_setfeasible( dim ) + set qh.FEASIBLEpoint from qh.feasible_string in "n,n,n" or "n n n" format + + notes: + "n,n,n" already checked by qh_initflags() + see qh_readfeasible() +*/ +void qh_setfeasible (int dim) { + int tokcount= 0; + char *s; + coordT *coords, value; + + if (!(s= qh feasible_string)) { + fprintf(qh ferr, "\ +qhull input error: halfspace intersection needs a feasible point.\n\ +Either prepend the input with 1 point or use 'Hn,n,n'. See manual.\n"); + qh_errexit (qh_ERRinput, NULL, NULL); + } + if (!(qh feasible_point= (pointT*)malloc (dim* sizeof(coordT)))) { + fprintf(qh ferr, "qhull error: insufficient memory for 'Hn,n,n'\n"); + qh_errexit(qh_ERRmem, NULL, NULL); + } + coords= qh feasible_point; + while (*s) { + value= qh_strtod (s, &s); + if (++tokcount > dim) { + fprintf (qh ferr, "qhull input warning: more coordinates for 'H%s' than dimension %d\n", + qh feasible_string, dim); + break; + } + *(coords++)= value; + if (*s) + s++; + } + while (++tokcount <= dim) + *(coords++)= 0.0; +} /* setfeasible */ + +/*-<a href="qh-io.htm#TOC" + >-------------------------------</a><a name="skipfacet">-</a> + + qh_skipfacet( facet ) + returns 'True' if this facet is not to be printed + + notes: + based on the user provided slice thresholds and 'good' specifications +*/ +boolT qh_skipfacet(facetT *facet) { + facetT *neighbor, **neighborp; + + if (qh PRINTneighbors) { + if (facet->good) + return !qh PRINTgood; + FOREACHneighbor_(facet) { + if (neighbor->good) + return False; + } + return True; + }else if (qh PRINTgood) + return !facet->good; + else if (!facet->normal) + return True; + return (!qh_inthresholds (facet->normal, NULL)); +} /* skipfacet */ + |