|
|
@@ -0,0 +1,4128 @@
|
|
|
+/*<html><pre> -<a href="qh-io_r.htm"
|
|
|
+ >-------------------------------</a><a name="TOP">-</a>
|
|
|
+
|
|
|
+ io_r.c
|
|
|
+ Input/Output routines of qhull application
|
|
|
+
|
|
|
+ see qh-io_r.htm and io_r.h
|
|
|
+
|
|
|
+ see user_r.c for qh_errprint and qh_printfacetlist
|
|
|
+
|
|
|
+ unix_r.c calls qh_readpoints and qh_produce_output
|
|
|
+
|
|
|
+ unix_r.c and user_r.c are the only callers of io_r.c functions
|
|
|
+ This allows the user to avoid loading io_r.o from qhull.a
|
|
|
+
|
|
|
+ Copyright (c) 1993-2020 The Geometry Center.
|
|
|
+ $Id: //main/2019/qhull/src/libqhull_r/io_r.c#12 $$Change: 2965 $
|
|
|
+ $DateTime: 2020/06/04 15:37:41 $$Author: bbarber $
|
|
|
+*/
|
|
|
+
|
|
|
+#include "qhull_ra.h"
|
|
|
+
|
|
|
+/*========= -functions in alphabetical order after qh_produce_output(qh) =====*/
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="produce_output">-</a>
|
|
|
+
|
|
|
+ qh_produce_output(qh )
|
|
|
+ qh_produce_output2(qh )
|
|
|
+ prints out the result of qhull in desired format
|
|
|
+ qh_produce_output2 does not call qh_prepare_output
|
|
|
+ qh_checkpolygon is valid for qh_prepare_output
|
|
|
+ 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(qhT *qh) {
|
|
|
+ int tempsize= qh_setsize(qh, qh->qhmem.tempstack);
|
|
|
+
|
|
|
+ qh_prepare_output(qh);
|
|
|
+ qh_produce_output2(qh);
|
|
|
+ if (qh_setsize(qh, qh->qhmem.tempstack) != tempsize) {
|
|
|
+ qh_fprintf(qh, qh->ferr, 6206, "qhull internal error (qh_produce_output): temporary sets not empty(%d)\n",
|
|
|
+ qh_setsize(qh, qh->qhmem.tempstack));
|
|
|
+ qh_errexit(qh, qh_ERRqhull, NULL, NULL);
|
|
|
+ }
|
|
|
+} /* produce_output */
|
|
|
+
|
|
|
+
|
|
|
+void qh_produce_output2(qhT *qh) {
|
|
|
+ int i, tempsize= qh_setsize(qh, qh->qhmem.tempstack), d_1;
|
|
|
+
|
|
|
+ fflush(NULL);
|
|
|
+ if (qh->PRINTsummary)
|
|
|
+ qh_printsummary(qh, qh->ferr);
|
|
|
+ else if (qh->PRINTout[0] == qh_PRINTnone)
|
|
|
+ qh_printsummary(qh, qh->fout);
|
|
|
+ for (i=0; i < qh_PRINTEND; i++)
|
|
|
+ qh_printfacets(qh, qh->fout, qh->PRINTout[i], qh->facet_list, NULL, !qh_ALL);
|
|
|
+ fflush(NULL);
|
|
|
+
|
|
|
+ qh_allstatistics(qh);
|
|
|
+ if (qh->PRINTprecision && !qh->MERGING && (qh->JOGGLEmax > REALmax/2 || qh->RERUN))
|
|
|
+ qh_printstats(qh, qh->ferr, qh->qhstat.precision, NULL);
|
|
|
+ if (qh->VERIFYoutput && (zzval_(Zridge) > 0 || zzval_(Zridgemid) > 0))
|
|
|
+ qh_printstats(qh, qh->ferr, qh->qhstat.vridges, NULL);
|
|
|
+ if (qh->PRINTstatistics) {
|
|
|
+ qh_printstatistics(qh, qh->ferr, "");
|
|
|
+ qh_memstatistics(qh, qh->ferr);
|
|
|
+ d_1= (int)sizeof(setT) + (qh->hull_dim - 1) * SETelemsize;
|
|
|
+ qh_fprintf(qh, qh->ferr, 8040, "\
|
|
|
+ size in bytes: merge %d ridge %d vertex %d facet %d\n\
|
|
|
+ normal %d ridge vertices %d facet vertices or neighbors %d\n",
|
|
|
+ (int)sizeof(mergeT), (int)sizeof(ridgeT),
|
|
|
+ (int)sizeof(vertexT), (int)sizeof(facetT),
|
|
|
+ qh->normal_size, d_1, d_1 + SETelemsize);
|
|
|
+ }
|
|
|
+ if (qh_setsize(qh, qh->qhmem.tempstack) != tempsize) {
|
|
|
+ qh_fprintf(qh, qh->ferr, 6065, "qhull internal error (qh_produce_output2): temporary sets not empty(%d)\n",
|
|
|
+ qh_setsize(qh, qh->qhmem.tempstack));
|
|
|
+ qh_errexit(qh, qh_ERRqhull, NULL, NULL);
|
|
|
+ }
|
|
|
+} /* produce_output2 */
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="dfacet">-</a>
|
|
|
+
|
|
|
+ qh_dfacet(qh, id )
|
|
|
+ print facet by id, for debugging
|
|
|
+
|
|
|
+*/
|
|
|
+void qh_dfacet(qhT *qh, unsigned int id) {
|
|
|
+ facetT *facet;
|
|
|
+
|
|
|
+ FORALLfacets {
|
|
|
+ if (facet->id == id) {
|
|
|
+ qh_printfacet(qh, qh->fout, facet);
|
|
|
+ break;
|
|
|
+ }
|
|
|
+ }
|
|
|
+} /* dfacet */
|
|
|
+
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="dvertex">-</a>
|
|
|
+
|
|
|
+ qh_dvertex(qh, id )
|
|
|
+ print vertex by id, for debugging
|
|
|
+*/
|
|
|
+void qh_dvertex(qhT *qh, unsigned int id) {
|
|
|
+ vertexT *vertex;
|
|
|
+
|
|
|
+ FORALLvertices {
|
|
|
+ if (vertex->id == id) {
|
|
|
+ qh_printvertex(qh, qh->fout, vertex);
|
|
|
+ break;
|
|
|
+ }
|
|
|
+ }
|
|
|
+} /* dvertex */
|
|
|
+
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="compare_facetarea">-</a>
|
|
|
+
|
|
|
+ qh_compare_facetarea( p1, p2 )
|
|
|
+ used by qsort() to order facets by area
|
|
|
+*/
|
|
|
+int qh_compare_facetarea(const void *p1, const void *p2) {
|
|
|
+ const facetT *a= *((facetT *const*)p1), *b= *((facetT *const*)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_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="compare_facetvisit">-</a>
|
|
|
+
|
|
|
+ qh_compare_facetvisit( p1, p2 )
|
|
|
+ used by qsort() to order facets by visit id or id
|
|
|
+*/
|
|
|
+int qh_compare_facetvisit(const void *p1, const void *p2) {
|
|
|
+ const facetT *a= *((facetT *const*)p1), *b= *((facetT *const*)p2);
|
|
|
+ int i,j;
|
|
|
+
|
|
|
+ if (!(i= (int)a->visitid))
|
|
|
+ i= (int)(0 - a->id); /* sign distinguishes id from visitid */
|
|
|
+ if (!(j= (int)b->visitid))
|
|
|
+ j= (int)(0 - b->id);
|
|
|
+ return(i - j);
|
|
|
+} /* compare_facetvisit */
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="compare_nummerge">-</a>
|
|
|
+
|
|
|
+ qh_compare_nummerge( p1, p2 )
|
|
|
+ used by qsort() to order facets by number of merges
|
|
|
+
|
|
|
+notes:
|
|
|
+ called by qh_markkeep ('PMerge-keep')
|
|
|
+*/
|
|
|
+int qh_compare_nummerge(const void *p1, const void *p2) {
|
|
|
+ const facetT *a= *((facetT *const*)p1), *b= *((facetT *const*)p2);
|
|
|
+
|
|
|
+ return(a->nummerge - b->nummerge);
|
|
|
+} /* compare_nummerge */
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="copyfilename">-</a>
|
|
|
+
|
|
|
+ qh_copyfilename(qh, dest, size, source, length )
|
|
|
+ copy filename identified by qh_skipfilename()
|
|
|
+
|
|
|
+ notes:
|
|
|
+ see qh_skipfilename() for syntax
|
|
|
+*/
|
|
|
+void qh_copyfilename(qhT *qh, char *filename, int size, const char* source, int length) {
|
|
|
+ char c= *source;
|
|
|
+
|
|
|
+ if (length > size + 1) {
|
|
|
+ qh_fprintf(qh, qh->ferr, 6040, "qhull error: filename is more than %d characters, %s\n", size-1, source);
|
|
|
+ qh_errexit(qh, qh_ERRinput, NULL, NULL);
|
|
|
+ }
|
|
|
+ strncpy(filename, source, (size_t)length);
|
|
|
+ filename[length]= '\0';
|
|
|
+ if (c == '\'' || c == '"') {
|
|
|
+ char *s= filename + 1;
|
|
|
+ char *t= filename;
|
|
|
+ while (*s) {
|
|
|
+ if (*s == c) {
|
|
|
+ if (s[-1] == '\\')
|
|
|
+ t[-1]= c;
|
|
|
+ }else
|
|
|
+ *t++= *s;
|
|
|
+ s++;
|
|
|
+ }
|
|
|
+ *t= '\0';
|
|
|
+ }
|
|
|
+} /* copyfilename */
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="countfacets">-</a>
|
|
|
+
|
|
|
+ qh_countfacets(qh, 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(qhT *qh, 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(qh, facet)))
|
|
|
+ facet->visitid= 0;
|
|
|
+ else {
|
|
|
+ facet->visitid= (unsigned int)(++numfacets);
|
|
|
+ totneighbors += qh_setsize(qh, facet->neighbors);
|
|
|
+ if (facet->simplicial) {
|
|
|
+ numsimplicial++;
|
|
|
+ if (facet->keepcentrum && facet->tricoplanar)
|
|
|
+ numtricoplanars++;
|
|
|
+ }else
|
|
|
+ numridges += qh_setsize(qh, facet->ridges);
|
|
|
+ if (facet->coplanarset)
|
|
|
+ numcoplanars += qh_setsize(qh, facet->coplanarset);
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ FOREACHfacet_(facets) {
|
|
|
+ if ((facet->visible && qh->NEWfacets)
|
|
|
+ || (!printall && qh_skipfacet(qh, facet)))
|
|
|
+ facet->visitid= 0;
|
|
|
+ else {
|
|
|
+ facet->visitid= (unsigned int)(++numfacets);
|
|
|
+ totneighbors += qh_setsize(qh, facet->neighbors);
|
|
|
+ if (facet->simplicial){
|
|
|
+ numsimplicial++;
|
|
|
+ if (facet->keepcentrum && facet->tricoplanar)
|
|
|
+ numtricoplanars++;
|
|
|
+ }else
|
|
|
+ numridges += qh_setsize(qh, facet->ridges);
|
|
|
+ if (facet->coplanarset)
|
|
|
+ numcoplanars += qh_setsize(qh, facet->coplanarset);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ qh->visit_id += (unsigned int)numfacets + 1;
|
|
|
+ *numfacetsp= numfacets;
|
|
|
+ *numsimplicialp= numsimplicial;
|
|
|
+ *totneighborsp= totneighbors;
|
|
|
+ *numridgesp= numridges;
|
|
|
+ *numcoplanarsp= numcoplanars;
|
|
|
+ *numtricoplanarsp= numtricoplanars;
|
|
|
+} /* countfacets */
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="detvnorm">-</a>
|
|
|
+
|
|
|
+ qh_detvnorm(qh, 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_r.c so that code for 'v Tv' can be removed by removing io_r.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(qhT *qh, 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, 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(qh, facet->vertices);
|
|
|
+ qh_setappend(qh, &points, facet->center);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ if (numcenters > dim) {
|
|
|
+ simplex= qh_settemp(qh, qh->TEMPsize);
|
|
|
+ qh_setappend(qh, &simplex, vertex->point);
|
|
|
+ if (unbounded)
|
|
|
+ qh_setappend(qh, &simplex, midpoint);
|
|
|
+ qh_maxsimplex(qh, dim, points, NULL, 0, &simplex);
|
|
|
+ qh_setdelnth(qh, simplex, 0);
|
|
|
+ }else if (numcenters == dim) {
|
|
|
+ if (unbounded)
|
|
|
+ qh_setappend(qh, &points, midpoint);
|
|
|
+ simplex= points;
|
|
|
+ }else {
|
|
|
+ qh_fprintf(qh, qh->ferr, 6216, "qhull internal error (qh_detvnorm): too few points(%d) to compute separating plane\n", numcenters);
|
|
|
+ qh_errexit(qh, qh_ERRqhull, NULL, NULL);
|
|
|
+ }
|
|
|
+ i= 0;
|
|
|
+ gmcoord= qh->gm_matrix;
|
|
|
+ point0= SETfirstt_(simplex, pointT);
|
|
|
+ FOREACHpoint_(simplex) {
|
|
|
+ if (qh->IStracing >= 4)
|
|
|
+ qh_printmatrix(qh, 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(qh, dim, qh->gm_row, point0, True,
|
|
|
+ normal, &offset, &nearzero);
|
|
|
+ /* nearzero is true for axis-parallel hyperplanes (e.g., a bounding box). Should detect degenerate hyperplanes. See 'Tv' check following */
|
|
|
+ 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(qh, vertex->point);
|
|
|
+ pointidA= qh_pointid(qh, 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, qh->ferr, 4014, "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(qh, 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, qh->ferr, 4015, "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_(qh, 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, qh->ferr, 4016, "qh_detvnorm: points %d %d Voronoi vertex %d dist %2.2g\n",
|
|
|
+ pointid, pointidA, facet->visitid, dist));
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ *offsetp= offset;
|
|
|
+ if (simplex != points)
|
|
|
+ qh_settempfree(qh, &simplex);
|
|
|
+ qh_settempfree(qh, &points);
|
|
|
+ return normal;
|
|
|
+} /* detvnorm */
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="detvridge">-</a>
|
|
|
+
|
|
|
+ qh_detvridge(qh, 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(qhT *qh, vertexT *vertex) {
|
|
|
+ setT *centers= qh_settemp(qh, qh->TEMPsize);
|
|
|
+ setT *tricenters= qh_settemp(qh, qh->TEMPsize);
|
|
|
+ facetT *neighbor, **neighborp;
|
|
|
+ boolT firstinf= True;
|
|
|
+
|
|
|
+ FOREACHneighbor_(vertex) {
|
|
|
+ if (neighbor->seen) {
|
|
|
+ if (neighbor->visitid) {
|
|
|
+ if (!neighbor->tricoplanar || qh_setunique(qh, &tricenters, neighbor->center))
|
|
|
+ qh_setappend(qh, ¢ers, neighbor);
|
|
|
+ }else if (firstinf) {
|
|
|
+ firstinf= False;
|
|
|
+ qh_setappend(qh, ¢ers, neighbor);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ qsort(SETaddr_(centers, facetT), (size_t)qh_setsize(qh, centers),
|
|
|
+ sizeof(facetT *), qh_compare_facetvisit);
|
|
|
+ qh_settempfree(qh, &tricenters);
|
|
|
+ return centers;
|
|
|
+} /* detvridge */
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="detvridge3">-</a>
|
|
|
+
|
|
|
+ qh_detvridge3(qh, 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 (!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(qhT *qh, vertexT *atvertex, vertexT *vertex) {
|
|
|
+ setT *centers= qh_settemp(qh, qh->TEMPsize);
|
|
|
+ setT *tricenters= qh_settemp(qh, 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(qh, &tricenters, facet->center))
|
|
|
+ qh_setappend(qh, ¢ers, facet);
|
|
|
+ }else if (firstinf) {
|
|
|
+ firstinf= False;
|
|
|
+ qh_setappend(qh, ¢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) {
|
|
|
+ qh_fprintf(qh, qh->ferr, 6217, "qhull internal error (qh_detvridge3): neighbors of vertex p%d are not connected at facet %d\n",
|
|
|
+ qh_pointid(qh, vertex->point), neighbor->id);
|
|
|
+ qh_errexit(qh, qh_ERRqhull, neighbor, NULL);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ FOREACHneighbor_(atvertex)
|
|
|
+ neighbor->seen2= True;
|
|
|
+ qh_settempfree(qh, &tricenters);
|
|
|
+ return centers;
|
|
|
+} /* detvridge3 */
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="eachvoronoi">-</a>
|
|
|
+
|
|
|
+ qh_eachvoronoi(qh, 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*)
|
|
|
+ fp may be NULL for QhullQh::qh_fprintf which calls appendQhullMessage
|
|
|
+ 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(qhT *qh, 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, 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(qh, tricenters, 0);
|
|
|
+ FOREACHneighborA_(vertex) {
|
|
|
+ if (neighborA->seen) {
|
|
|
+ if (neighborA->visitid) {
|
|
|
+ if (!neighborA->tricoplanar || qh_setunique(qh, &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, qh->ferr, 4017, "qh_eachvoronoi: Voronoi ridge of %d vertices between sites %d and %d\n",
|
|
|
+ count, qh_pointid(qh, atvertex->point), qh_pointid(qh, vertex->point)));
|
|
|
+ if (printvridge) {
|
|
|
+ if (inorder && qh->hull_dim == 3+1) /* 3-d Voronoi diagram */
|
|
|
+ centers= qh_detvridge3(qh, atvertex, vertex);
|
|
|
+ else
|
|
|
+ centers= qh_detvridge(qh, vertex);
|
|
|
+ (*printvridge)(qh, fp, atvertex, vertex, centers, unbounded);
|
|
|
+ qh_settempfree(qh, ¢ers);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ FOREACHneighbor_(atvertex)
|
|
|
+ neighbor->seen= False;
|
|
|
+ qh_settempfree(qh, &tricenters);
|
|
|
+ return totridges;
|
|
|
+} /* eachvoronoi */
|
|
|
+
|
|
|
+
|
|
|
+/*-<a href="qh-poly_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="eachvoronoi_all">-</a>
|
|
|
+
|
|
|
+ qh_eachvoronoi_all(qh, 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(qhT *qh, 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, qh_ASvoronoi);
|
|
|
+ qh_vertexneighbors(qh);
|
|
|
+ maximize_(qh->visit_id, (unsigned int)qh->num_facets);
|
|
|
+ FORALLfacets {
|
|
|
+ facet->visitid= 0;
|
|
|
+ facet->seen= False;
|
|
|
+ facet->seen2= True;
|
|
|
+ }
|
|
|
+ FORALLfacets {
|
|
|
+ if (facet->upperdelaunay == isUpper)
|
|
|
+ facet->visitid= (unsigned int)(numcenters++);
|
|
|
+ }
|
|
|
+ FORALLvertices
|
|
|
+ vertex->seen= False;
|
|
|
+ FORALLvertices {
|
|
|
+ if (qh->GOODvertex > 0 && qh_pointid(qh, vertex->point)+1 != qh->GOODvertex)
|
|
|
+ continue;
|
|
|
+ totridges += qh_eachvoronoi(qh, fp, printvridge, vertex,
|
|
|
+ !qh_ALL, innerouter, inorder);
|
|
|
+ }
|
|
|
+ return totridges;
|
|
|
+} /* eachvoronoi_all */
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="facet2point">-</a>
|
|
|
+
|
|
|
+ qh_facet2point(qh, 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(qhT *qh, 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(qh, vertex0->point, facet, &dist);
|
|
|
+ *mindist= dist;
|
|
|
+ *point0= qh_projectpoint(qh, vertex0->point, facet, dist);
|
|
|
+ qh_distplane(qh, vertex1->point, facet, &dist);
|
|
|
+ minimize_(*mindist, dist);
|
|
|
+ *point1= qh_projectpoint(qh, vertex1->point, facet, dist);
|
|
|
+} /* facet2point */
|
|
|
+
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="facetvertices">-</a>
|
|
|
+
|
|
|
+ qh_facetvertices(qh, 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(qhT *qh, 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, qh->num_vertices);
|
|
|
+ FORALLvertices {
|
|
|
+ vertex->visitid= qh->vertex_visit;
|
|
|
+ qh_setappend(qh, &vertices, vertex);
|
|
|
+ }
|
|
|
+ }else {
|
|
|
+ vertices= qh_settemp(qh, qh->TEMPsize);
|
|
|
+ FORALLfacet_(facetlist) {
|
|
|
+ if (!allfacets && qh_skipfacet(qh, facet))
|
|
|
+ continue;
|
|
|
+ FOREACHvertex_(facet->vertices) {
|
|
|
+ if (vertex->visitid != qh->vertex_visit) {
|
|
|
+ vertex->visitid= qh->vertex_visit;
|
|
|
+ qh_setappend(qh, &vertices, vertex);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ FOREACHfacet_(facets) {
|
|
|
+ if (!allfacets && qh_skipfacet(qh, facet))
|
|
|
+ continue;
|
|
|
+ FOREACHvertex_(facet->vertices) {
|
|
|
+ if (vertex->visitid != qh->vertex_visit) {
|
|
|
+ vertex->visitid= qh->vertex_visit;
|
|
|
+ qh_setappend(qh, &vertices, vertex);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ return vertices;
|
|
|
+} /* facetvertices */
|
|
|
+
|
|
|
+/*-<a href="qh-geom_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="geomplanes">-</a>
|
|
|
+
|
|
|
+ qh_geomplanes(qh, 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_r.c with roundoff covered by qh_GEOMepsilon
|
|
|
+*/
|
|
|
+void qh_geomplanes(qhT *qh, facetT *facet, realT *outerplane, realT *innerplane) {
|
|
|
+ realT radius;
|
|
|
+
|
|
|
+ if (qh->MERGING || qh->JOGGLEmax < REALmax/2) {
|
|
|
+ qh_outerinner(qh, facet, outerplane, innerplane);
|
|
|
+ radius= qh->PRINTradius;
|
|
|
+ if (qh->JOGGLEmax < REALmax/2)
|
|
|
+ radius -= qh->JOGGLEmax * sqrt((realT)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_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="markkeep">-</a>
|
|
|
+
|
|
|
+ qh_markkeep(qh, facetlist )
|
|
|
+ restrict good facets for qh.KEEParea, qh.KEEPmerge, and qh.KEEPminArea
|
|
|
+ ignores visible facets (!part of convex hull)
|
|
|
+
|
|
|
+ returns:
|
|
|
+ may clear facet->good
|
|
|
+ recomputes qh.num_good
|
|
|
+
|
|
|
+ notes:
|
|
|
+ only called by qh_prepare_output after qh_findgood_all
|
|
|
+ does not throw errors except memory/corruption of qset_r.c
|
|
|
+
|
|
|
+ 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(qhT *qh, facetT *facetlist) {
|
|
|
+ facetT *facet, **facetp;
|
|
|
+ setT *facets= qh_settemp(qh, qh->num_facets);
|
|
|
+ int size, count;
|
|
|
+
|
|
|
+ trace2((qh, qh->ferr, 2006, "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(qh, &facets, facet);
|
|
|
+ }
|
|
|
+ size= qh_setsize(qh, facets);
|
|
|
+ if (qh->KEEParea) {
|
|
|
+ qsort(SETaddr_(facets, facetT), (size_t)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_t)size,
|
|
|
+ sizeof(facetT *), qh_compare_nummerge);
|
|
|
+ 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(qh, &facets);
|
|
|
+ count= 0;
|
|
|
+ FORALLfacet_(facetlist) {
|
|
|
+ if (facet->good)
|
|
|
+ count++;
|
|
|
+ }
|
|
|
+ qh->num_good= count;
|
|
|
+} /* markkeep */
|
|
|
+
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="markvoronoi">-</a>
|
|
|
+
|
|
|
+ qh_markvoronoi(qh, 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(qhT *qh, 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, qh_ASvoronoi); /* in case, qh_printvdiagram2 called by user */
|
|
|
+ qh_vertexneighbors(qh);
|
|
|
+ vertices= qh_pointvertex(qh);
|
|
|
+ if (qh->ATinfinity)
|
|
|
+ SETelem_(vertices, qh->num_points-1)= NULL;
|
|
|
+ qh->visit_id++;
|
|
|
+ maximize_(qh->visit_id, (unsigned int)qh->num_facets);
|
|
|
+ FORALLfacet_(facetlist) {
|
|
|
+ if (printall || !qh_skipfacet(qh, facet)) {
|
|
|
+ if (!facet->upperdelaunay) {
|
|
|
+ isLower= True;
|
|
|
+ break;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ FOREACHfacet_(facets) {
|
|
|
+ if (printall || !qh_skipfacet(qh, 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(qh, facet))
|
|
|
+ facet->visitid= (unsigned int)(numcenters++);
|
|
|
+ }
|
|
|
+ FOREACHfacet_(facets) {
|
|
|
+ if (printall || !qh_skipfacet(qh, facet))
|
|
|
+ facet->visitid= (unsigned int)(numcenters++);
|
|
|
+ }
|
|
|
+ *isLowerp= isLower;
|
|
|
+ *numcentersp= numcenters;
|
|
|
+ trace2((qh, qh->ferr, 2007, "qh_markvoronoi: isLower %d numcenters %d\n", isLower, numcenters));
|
|
|
+ return vertices;
|
|
|
+} /* markvoronoi */
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="order_vertexneighbors">-</a>
|
|
|
+
|
|
|
+ qh_order_vertexneighbors(qh, vertex )
|
|
|
+ order facet neighbors of vertex by 2-d (orientation), 3-d (adjacency), or n-d (f.visitid,id)
|
|
|
+
|
|
|
+ notes:
|
|
|
+ error if qh_vertexneighbors not called beforehand
|
|
|
+ only 2-d orients the neighbors
|
|
|
+ for 4-d and higher
|
|
|
+ set or clear f.visitid for qh_compare_facetvisit
|
|
|
+ for example, use qh_markvoronoi (e.g., qh_printvornoi) or qh_countfacets (e.g., qh_printvneighbors)
|
|
|
+
|
|
|
+ design (2-d):
|
|
|
+ see qh_printextremes_2d
|
|
|
+ design (3-d):
|
|
|
+ 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
|
|
|
+ design (n-d):
|
|
|
+ qsort by f.visitid, or f.facetid (qh_compare_facetvisit)
|
|
|
+ facet_id is negated (sorted before visit_id facets)
|
|
|
+*/
|
|
|
+void qh_order_vertexneighbors(qhT *qh, vertexT *vertex) {
|
|
|
+ setT *newset;
|
|
|
+ facetT *facet, *facetA, *facetB, *neighbor, **neighborp;
|
|
|
+ vertexT *vertexA;
|
|
|
+ int numneighbors;
|
|
|
+
|
|
|
+ trace4((qh, qh->ferr, 4018, "qh_order_vertexneighbors: order facet neighbors of v%d by 2-d (orientation), 3-d (adjacency), or n-d (f.visitid,id)\n", vertex->id));
|
|
|
+ if (!qh->VERTEXneighbors) {
|
|
|
+ qh_fprintf(qh, qh->ferr, 6428, "qhull internal error (qh_order_vertexneighbors): call qh_vertexneighbors before calling qh_order_vertexneighbors\n");
|
|
|
+ qh_errexit(qh, qh_ERRqhull, NULL, NULL);
|
|
|
+ }
|
|
|
+ if (qh->hull_dim == 2) {
|
|
|
+ facetA= SETfirstt_(vertex->neighbors, facetT);
|
|
|
+ if (facetA->toporient ^ qh_ORIENTclock)
|
|
|
+ vertexA= SETfirstt_(facetA->vertices, vertexT);
|
|
|
+ else
|
|
|
+ vertexA= SETsecondt_(facetA->vertices, vertexT);
|
|
|
+ if (vertexA!=vertex) {
|
|
|
+ facetB= SETsecondt_(vertex->neighbors, facetT);
|
|
|
+ SETfirst_(vertex->neighbors)= facetB;
|
|
|
+ SETsecond_(vertex->neighbors)= facetA;
|
|
|
+ }
|
|
|
+ }else if (qh->hull_dim == 3) {
|
|
|
+ newset= qh_settemp(qh, qh_setsize(qh, vertex->neighbors));
|
|
|
+ facet= (facetT *)qh_setdellast(vertex->neighbors);
|
|
|
+ qh_setappend(qh, &newset, facet);
|
|
|
+ while (qh_setsize(qh, vertex->neighbors)) {
|
|
|
+ FOREACHneighbor_(vertex) {
|
|
|
+ if (qh_setin(facet->neighbors, neighbor)) {
|
|
|
+ qh_setdel(vertex->neighbors, neighbor);
|
|
|
+ qh_setappend(qh, &newset, neighbor);
|
|
|
+ facet= neighbor;
|
|
|
+ break;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ if (!neighbor) {
|
|
|
+ qh_fprintf(qh, qh->ferr, 6066, "qhull internal error (qh_order_vertexneighbors): no neighbor of v%d for f%d\n",
|
|
|
+ vertex->id, facet->id);
|
|
|
+ qh_errexit(qh, qh_ERRqhull, facet, NULL);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ qh_setfree(qh, &vertex->neighbors);
|
|
|
+ qh_settemppop(qh);
|
|
|
+ vertex->neighbors= newset;
|
|
|
+ }else { /* qh.hull_dim >= 4 */
|
|
|
+ numneighbors= qh_setsize(qh, vertex->neighbors);
|
|
|
+ qsort(SETaddr_(vertex->neighbors, facetT), (size_t)numneighbors,
|
|
|
+ sizeof(facetT *), qh_compare_facetvisit);
|
|
|
+ }
|
|
|
+} /* order_vertexneighbors */
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="prepare_output">-</a>
|
|
|
+
|
|
|
+ qh_prepare_output(qh )
|
|
|
+ prepare for qh_produce_output2(qh) according to
|
|
|
+ qh.KEEPminArea, KEEParea, KEEPmerge, GOODvertex, GOODthreshold, GOODpoint, ONLYgood, SPLITthresholds
|
|
|
+ does not reset facet->good
|
|
|
+
|
|
|
+ notes
|
|
|
+ called by qh_produce_output, qh_new_qhull, Qhull.outputQhull
|
|
|
+ except for PRINTstatistics, no-op if previously called with same options
|
|
|
+*/
|
|
|
+void qh_prepare_output(qhT *qh) {
|
|
|
+ if (qh->VORONOI) {
|
|
|
+ qh_clearcenters(qh, qh_ASvoronoi); /* must be before qh_triangulate */
|
|
|
+ qh_vertexneighbors(qh);
|
|
|
+ }
|
|
|
+ if (qh->TRIangulate && !qh->hasTriangulation) {
|
|
|
+ qh_triangulate(qh);
|
|
|
+ if (qh->VERIFYoutput && !qh->CHECKfrequently)
|
|
|
+ qh_checkpolygon(qh, qh->facet_list);
|
|
|
+ }
|
|
|
+ qh_findgood_all(qh, qh->facet_list);
|
|
|
+ if (qh->GETarea)
|
|
|
+ qh_getarea(qh, qh->facet_list);
|
|
|
+ if (qh->KEEParea || qh->KEEPmerge || qh->KEEPminArea < REALmax/2)
|
|
|
+ qh_markkeep(qh, qh->facet_list);
|
|
|
+ if (qh->PRINTstatistics)
|
|
|
+ qh_collectstatistics(qh);
|
|
|
+}
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printafacet">-</a>
|
|
|
+
|
|
|
+ qh_printafacet(qh, 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(qhT *qh, FILE *fp, qh_PRINT 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(qh, facet))
|
|
|
+ return;
|
|
|
+ if (facet->visible && qh->NEWfacets && format != qh_PRINTfacets)
|
|
|
+ return;
|
|
|
+ qh->printoutnum++;
|
|
|
+ switch (format) {
|
|
|
+ case qh_PRINTarea:
|
|
|
+ if (facet->isarea) {
|
|
|
+ qh_fprintf(qh, fp, 9009, qh_REAL_1, facet->f.area);
|
|
|
+ qh_fprintf(qh, fp, 9010, "\n");
|
|
|
+ }else
|
|
|
+ qh_fprintf(qh, fp, 9011, "0\n");
|
|
|
+ break;
|
|
|
+ case qh_PRINTcoplanars:
|
|
|
+ qh_fprintf(qh, fp, 9012, "%d", qh_setsize(qh, facet->coplanarset));
|
|
|
+ FOREACHpoint_(facet->coplanarset)
|
|
|
+ qh_fprintf(qh, fp, 9013, " %d", qh_pointid(qh, point));
|
|
|
+ qh_fprintf(qh, fp, 9014, "\n");
|
|
|
+ break;
|
|
|
+ case qh_PRINTcentrums:
|
|
|
+ qh_printcenter(qh, fp, format, NULL, facet);
|
|
|
+ break;
|
|
|
+ case qh_PRINTfacets:
|
|
|
+ qh_printfacet(qh, fp, facet);
|
|
|
+ break;
|
|
|
+ case qh_PRINTfacets_xridge:
|
|
|
+ qh_printfacetheader(qh, 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(qh, color, color);
|
|
|
+ if (qh->PRINTdim != qh->hull_dim)
|
|
|
+ qh_normalize2(qh, color, 3, True, NULL, NULL);
|
|
|
+ if (qh->hull_dim <= 2)
|
|
|
+ qh_printfacet2geom(qh, fp, facet, color);
|
|
|
+ else if (qh->hull_dim == 3) {
|
|
|
+ if (facet->simplicial)
|
|
|
+ qh_printfacet3geom_simplicial(qh, fp, facet, color);
|
|
|
+ else
|
|
|
+ qh_printfacet3geom_nonsimplicial(qh, fp, facet, color);
|
|
|
+ }else {
|
|
|
+ if (facet->simplicial)
|
|
|
+ qh_printfacet4geom_simplicial(qh, fp, facet, color);
|
|
|
+ else
|
|
|
+ qh_printfacet4geom_nonsimplicial(qh, fp, facet, color);
|
|
|
+ }
|
|
|
+ break;
|
|
|
+ case qh_PRINTids:
|
|
|
+ qh_fprintf(qh, fp, 9015, "%d\n", facet->id);
|
|
|
+ break;
|
|
|
+ case qh_PRINTincidences:
|
|
|
+ case qh_PRINToff:
|
|
|
+ case qh_PRINTtriangles:
|
|
|
+ if (qh->hull_dim == 3 && format != qh_PRINTtriangles)
|
|
|
+ qh_printfacet3vertex(qh, fp, facet, format);
|
|
|
+ else if (facet->simplicial || qh->hull_dim == 2 || format == qh_PRINToff)
|
|
|
+ qh_printfacetNvertex_simplicial(qh, fp, facet, format);
|
|
|
+ else
|
|
|
+ qh_printfacetNvertex_nonsimplicial(qh, fp, facet, qh->printoutvar++, format);
|
|
|
+ break;
|
|
|
+ case qh_PRINTinner:
|
|
|
+ qh_outerinner(qh, facet, NULL, &innerplane);
|
|
|
+ offset= facet->offset - innerplane;
|
|
|
+ goto LABELprintnorm;
|
|
|
+ break; /* prevent warning */
|
|
|
+ case qh_PRINTmerges:
|
|
|
+ qh_fprintf(qh, fp, 9016, "%d\n", facet->nummerge);
|
|
|
+ break;
|
|
|
+ case qh_PRINTnormals:
|
|
|
+ offset= facet->offset;
|
|
|
+ goto LABELprintnorm;
|
|
|
+ break; /* prevent warning */
|
|
|
+ case qh_PRINTouter:
|
|
|
+ qh_outerinner(qh, facet, &outerplane, NULL);
|
|
|
+ offset= facet->offset - outerplane;
|
|
|
+ LABELprintnorm:
|
|
|
+ if (!facet->normal) {
|
|
|
+ qh_fprintf(qh, fp, 9017, "no normal for facet f%d\n", facet->id);
|
|
|
+ break;
|
|
|
+ }
|
|
|
+ if (qh->CDDoutput) {
|
|
|
+ qh_fprintf(qh, fp, 9018, qh_REAL_1, -offset);
|
|
|
+ for (k=0; k < qh->hull_dim; k++)
|
|
|
+ qh_fprintf(qh, fp, 9019, qh_REAL_1, -facet->normal[k]);
|
|
|
+ }else {
|
|
|
+ for (k=0; k < qh->hull_dim; k++)
|
|
|
+ qh_fprintf(qh, fp, 9020, qh_REAL_1, facet->normal[k]);
|
|
|
+ qh_fprintf(qh, fp, 9021, qh_REAL_1, offset);
|
|
|
+ }
|
|
|
+ qh_fprintf(qh, fp, 9022, "\n");
|
|
|
+ break;
|
|
|
+ case qh_PRINTmathematica: /* either 2 or 3-d by qh_printbegin */
|
|
|
+ case qh_PRINTmaple:
|
|
|
+ if (qh->hull_dim == 2)
|
|
|
+ qh_printfacet2math(qh, fp, facet, format, qh->printoutvar++);
|
|
|
+ else
|
|
|
+ qh_printfacet3math(qh, fp, facet, format, qh->printoutvar++);
|
|
|
+ break;
|
|
|
+ case qh_PRINTneighbors:
|
|
|
+ qh_fprintf(qh, fp, 9023, "%d", qh_setsize(qh, facet->neighbors));
|
|
|
+ FOREACHneighbor_(facet)
|
|
|
+ qh_fprintf(qh, fp, 9024, " %d",
|
|
|
+ neighbor->visitid ? neighbor->visitid - 1: 0 - neighbor->id);
|
|
|
+ qh_fprintf(qh, fp, 9025, "\n");
|
|
|
+ break;
|
|
|
+ case qh_PRINTpointintersect:
|
|
|
+ if (!qh->feasible_point) {
|
|
|
+ qh_fprintf(qh, qh->ferr, 6067, "qhull input error (qh_printafacet): option 'Fp' needs qh->feasible_point\n");
|
|
|
+ qh_errexit(qh, qh_ERRinput, NULL, NULL);
|
|
|
+ }
|
|
|
+ if (facet->offset > 0)
|
|
|
+ goto LABELprintinfinite;
|
|
|
+ point= coordp= (coordT *)qh_memalloc(qh, 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(qh, point, qh->normal_size);
|
|
|
+ goto LABELprintinfinite;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ qh_printpoint(qh, fp, NULL, point);
|
|
|
+ qh_memfree(qh, point, qh->normal_size);
|
|
|
+ break;
|
|
|
+ LABELprintinfinite:
|
|
|
+ for (k=qh->hull_dim; k--; )
|
|
|
+ qh_fprintf(qh, fp, 9026, qh_REAL_1, qh_INFINITE);
|
|
|
+ qh_fprintf(qh, fp, 9027, "\n");
|
|
|
+ break;
|
|
|
+ case qh_PRINTpointnearest:
|
|
|
+ FOREACHpoint_(facet->coplanarset) {
|
|
|
+ int id, id2;
|
|
|
+ vertex= qh_nearvertex(qh, facet, point, &dist);
|
|
|
+ id= qh_pointid(qh, vertex->point);
|
|
|
+ id2= qh_pointid(qh, point);
|
|
|
+ qh_fprintf(qh, fp, 9028, "%d %d %d " qh_REAL_1 "\n", id, id2, facet->id, dist);
|
|
|
+ }
|
|
|
+ break;
|
|
|
+ case qh_PRINTpoints: /* VORONOI only by qh_printbegin */
|
|
|
+ if (qh->CDDoutput)
|
|
|
+ qh_fprintf(qh, fp, 9029, "1 ");
|
|
|
+ qh_printcenter(qh, fp, format, NULL, facet);
|
|
|
+ break;
|
|
|
+ case qh_PRINTvertices:
|
|
|
+ qh_fprintf(qh, fp, 9030, "%d", qh_setsize(qh, facet->vertices));
|
|
|
+ FOREACHvertex_(facet->vertices)
|
|
|
+ qh_fprintf(qh, fp, 9031, " %d", qh_pointid(qh, vertex->point));
|
|
|
+ qh_fprintf(qh, fp, 9032, "\n");
|
|
|
+ break;
|
|
|
+ default:
|
|
|
+ break;
|
|
|
+ }
|
|
|
+} /* printafacet */
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printbegin">-</a>
|
|
|
+
|
|
|
+ qh_printbegin(qh )
|
|
|
+ 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(qhT *qh, FILE *fp, qh_PRINT 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(qh, facetlist, facets, printall, &numfacets, &numsimplicial,
|
|
|
+ &totneighbors, &numridges, &numcoplanars, &numtricoplanars);
|
|
|
+ switch (format) {
|
|
|
+ case qh_PRINTnone:
|
|
|
+ break;
|
|
|
+ case qh_PRINTarea:
|
|
|
+ qh_fprintf(qh, fp, 9033, "%d\n", numfacets);
|
|
|
+ break;
|
|
|
+ case qh_PRINTcoplanars:
|
|
|
+ qh_fprintf(qh, fp, 9034, "%d\n", numfacets);
|
|
|
+ break;
|
|
|
+ case qh_PRINTcentrums:
|
|
|
+ if (qh->CENTERtype == qh_ASnone)
|
|
|
+ qh_clearcenters(qh, qh_AScentrum);
|
|
|
+ qh_fprintf(qh, fp, 9035, "%d\n%d\n", qh->hull_dim, numfacets);
|
|
|
+ break;
|
|
|
+ case qh_PRINTfacets:
|
|
|
+ case qh_PRINTfacets_xridge:
|
|
|
+ if (facetlist)
|
|
|
+ qh_printvertexlist(qh, 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))
|
|
|
+ qh_fprintf(qh, qh->ferr, 7049, "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)))
|
|
|
+ qh_fprintf(qh, qh->ferr, 7050, "qhull warning: output for outer/inner planes and centrums not implemented in 4-d\n");
|
|
|
+ if (qh->PRINTdim == 4 && (qh->PRINTspheres))
|
|
|
+ qh_fprintf(qh, qh->ferr, 7051, "qhull warning: output for vertices not implemented in 4-d\n");
|
|
|
+ if (qh->PRINTdim == 4 && qh->DOintersections && qh->PRINTnoplanes)
|
|
|
+ qh_fprintf(qh, qh->ferr, 7052, "qhull warning: 'Gnh' generates no output in 4-d\n");
|
|
|
+ if (qh->PRINTdim == 2) {
|
|
|
+ qh_fprintf(qh, fp, 9036, "{appearance {linewidth 3} LIST # %s | %s\n",
|
|
|
+ qh->rbox_command, qh->qhull_command);
|
|
|
+ }else if (qh->PRINTdim == 3) {
|
|
|
+ qh_fprintf(qh, fp, 9037, "{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(qh, NULL, facet, &num, printall);
|
|
|
+ FOREACHfacet_(facets)
|
|
|
+ qh_printend4geom(qh, NULL, facet, &num, printall);
|
|
|
+ qh->ridgeoutnum= num;
|
|
|
+ qh->printoutvar= 0; /* counts number of ridges in output */
|
|
|
+ qh_fprintf(qh, fp, 9038, "LIST # %s | %s\n", qh->rbox_command, qh->qhull_command);
|
|
|
+ }
|
|
|
+
|
|
|
+ if (qh->PRINTdots) {
|
|
|
+ qh->printoutnum++;
|
|
|
+ num= qh->num_points + qh_setsize(qh, qh->other_points);
|
|
|
+ if (qh->DELAUNAY && qh->ATinfinity)
|
|
|
+ num--;
|
|
|
+ if (qh->PRINTdim == 4)
|
|
|
+ qh_fprintf(qh, fp, 9039, "4VECT %d %d 1\n", num, num);
|
|
|
+ else
|
|
|
+ qh_fprintf(qh, fp, 9040, "VECT %d %d 1\n", num, num);
|
|
|
+
|
|
|
+ for (i=num; i--; ) {
|
|
|
+ if (i % 20 == 0)
|
|
|
+ qh_fprintf(qh, fp, 9041, "\n");
|
|
|
+ qh_fprintf(qh, fp, 9042, "1 ");
|
|
|
+ }
|
|
|
+ qh_fprintf(qh, fp, 9043, "# 1 point per line\n1 ");
|
|
|
+ for (i=num-1; i--; ) { /* num at least 3 for D2 */
|
|
|
+ if (i % 20 == 0)
|
|
|
+ qh_fprintf(qh, fp, 9044, "\n");
|
|
|
+ qh_fprintf(qh, fp, 9045, "0 ");
|
|
|
+ }
|
|
|
+ qh_fprintf(qh, fp, 9046, "# 1 color for all\n");
|
|
|
+ FORALLpoints {
|
|
|
+ if (!qh->DELAUNAY || !qh->ATinfinity || qh_pointid(qh, point) != qh->num_points-1) {
|
|
|
+ if (qh->PRINTdim == 4)
|
|
|
+ qh_printpoint(qh, fp, NULL, point);
|
|
|
+ else
|
|
|
+ qh_printpoint3(qh, fp, point);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ FOREACHpoint_(qh->other_points) {
|
|
|
+ if (qh->PRINTdim == 4)
|
|
|
+ qh_printpoint(qh, fp, NULL, point);
|
|
|
+ else
|
|
|
+ qh_printpoint3(qh, fp, point);
|
|
|
+ }
|
|
|
+ qh_fprintf(qh, fp, 9047, "0 1 1 1 # color of points\n");
|
|
|
+ }
|
|
|
+
|
|
|
+ if (qh->PRINTdim == 4 && !qh->PRINTnoplanes)
|
|
|
+ /* 4dview loads up multiple 4OFF objects slowly */
|
|
|
+ qh_fprintf(qh, fp, 9048, "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((realT)qh->hull_dim);
|
|
|
+ if (qh->PRINTdim != 4 &&
|
|
|
+ (qh->PRINTcoplanar || qh->PRINTspheres || qh->PRINTcentrums)) {
|
|
|
+ vertices= qh_facetvertices(qh, facetlist, facets, printall);
|
|
|
+ if (qh->PRINTspheres && qh->PRINTdim <= 3)
|
|
|
+ qh_printspheres(qh, fp, vertices, qh->PRINTradius);
|
|
|
+ if (qh->PRINTcoplanar || qh->PRINTcentrums) {
|
|
|
+ qh->firstcentrum= True;
|
|
|
+ if (qh->PRINTcoplanar&& !qh->PRINTspheres) {
|
|
|
+ FOREACHvertex_(vertices)
|
|
|
+ qh_printpointvect2(qh, fp, vertex->point, NULL, qh->interior_point, qh->PRINTradius);
|
|
|
+ }
|
|
|
+ FORALLfacet_(facetlist) {
|
|
|
+ if (!printall && qh_skipfacet(qh, facet))
|
|
|
+ continue;
|
|
|
+ if (!facet->normal)
|
|
|
+ continue;
|
|
|
+ if (qh->PRINTcentrums && qh->PRINTdim <= 3)
|
|
|
+ qh_printcentrum(qh, fp, facet, qh->PRINTcradius);
|
|
|
+ if (!qh->PRINTcoplanar)
|
|
|
+ continue;
|
|
|
+ FOREACHpoint_(facet->coplanarset)
|
|
|
+ qh_printpointvect2(qh, fp, point, facet->normal, NULL, qh->PRINTradius);
|
|
|
+ FOREACHpoint_(facet->outsideset)
|
|
|
+ qh_printpointvect2(qh, fp, point, facet->normal, NULL, qh->PRINTradius);
|
|
|
+ }
|
|
|
+ FOREACHfacet_(facets) {
|
|
|
+ if (!printall && qh_skipfacet(qh, facet))
|
|
|
+ continue;
|
|
|
+ if (!facet->normal)
|
|
|
+ continue;
|
|
|
+ if (qh->PRINTcentrums && qh->PRINTdim <= 3)
|
|
|
+ qh_printcentrum(qh, fp, facet, qh->PRINTcradius);
|
|
|
+ if (!qh->PRINTcoplanar)
|
|
|
+ continue;
|
|
|
+ FOREACHpoint_(facet->coplanarset)
|
|
|
+ qh_printpointvect2(qh, fp, point, facet->normal, NULL, qh->PRINTradius);
|
|
|
+ FOREACHpoint_(facet->outsideset)
|
|
|
+ qh_printpointvect2(qh, fp, point, facet->normal, NULL, qh->PRINTradius);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ qh_settempfree(qh, &vertices);
|
|
|
+ }
|
|
|
+ qh->visit_id++; /* for printing hyperplane intersections */
|
|
|
+ break;
|
|
|
+ case qh_PRINTids:
|
|
|
+ qh_fprintf(qh, fp, 9049, "%d\n", numfacets);
|
|
|
+ break;
|
|
|
+ case qh_PRINTincidences:
|
|
|
+ if (qh->VORONOI && qh->PRINTprecision)
|
|
|
+ qh_fprintf(qh, qh->ferr, 7053, "qhull warning: input sites of Delaunay regions (option 'i'). Use option 'p' or 'o' for Voronoi centers. Disable warning with option 'Pp'\n");
|
|
|
+ qh->printoutvar= (int)qh->vertex_id; /* centrum id for 4-d+, non-simplicial facets */
|
|
|
+ if (qh->hull_dim <= 3)
|
|
|
+ qh_fprintf(qh, fp, 9050, "%d\n", numfacets);
|
|
|
+ else
|
|
|
+ qh_fprintf(qh, fp, 9051, "%d\n", numsimplicial+numridges);
|
|
|
+ break;
|
|
|
+ case qh_PRINTinner:
|
|
|
+ case qh_PRINTnormals:
|
|
|
+ case qh_PRINTouter:
|
|
|
+ if (qh->CDDoutput)
|
|
|
+ qh_fprintf(qh, fp, 9052, "%s | %s\nbegin\n %d %d real\n", qh->rbox_command,
|
|
|
+ qh->qhull_command, numfacets, qh->hull_dim+1);
|
|
|
+ else
|
|
|
+ qh_fprintf(qh, fp, 9053, "%d\n%d\n", qh->hull_dim+1, numfacets);
|
|
|
+ break;
|
|
|
+ case qh_PRINTmathematica:
|
|
|
+ case qh_PRINTmaple:
|
|
|
+ if (qh->hull_dim > 3) /* qh_initbuffers also checks */
|
|
|
+ goto LABELnoformat;
|
|
|
+ if (qh->VORONOI)
|
|
|
+ qh_fprintf(qh, qh->ferr, 7054, "qhull warning: output is the Delaunay triangulation\n");
|
|
|
+ if (format == qh_PRINTmaple) {
|
|
|
+ if (qh->hull_dim == 2)
|
|
|
+ qh_fprintf(qh, fp, 9054, "PLOT(CURVES(\n");
|
|
|
+ else
|
|
|
+ qh_fprintf(qh, fp, 9055, "PLOT3D(POLYGONS(\n");
|
|
|
+ }else
|
|
|
+ qh_fprintf(qh, fp, 9056, "{\n");
|
|
|
+ qh->printoutvar= 0; /* counts number of facets for notfirst */
|
|
|
+ break;
|
|
|
+ case qh_PRINTmerges:
|
|
|
+ qh_fprintf(qh, fp, 9057, "%d\n", numfacets);
|
|
|
+ break;
|
|
|
+ case qh_PRINTpointintersect:
|
|
|
+ qh_fprintf(qh, fp, 9058, "%d\n%d\n", qh->hull_dim, numfacets);
|
|
|
+ break;
|
|
|
+ case qh_PRINTneighbors:
|
|
|
+ qh_fprintf(qh, fp, 9059, "%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)
|
|
|
+ qh_fprintf(qh, fp, 9060, "%d\n%d %d %d\n", num,
|
|
|
+ qh->num_points+qh_setsize(qh, qh->other_points), numfacets, totneighbors/2);
|
|
|
+ else { /* qh_PRINTtriangles */
|
|
|
+ qh->printoutvar= qh->num_points+qh_setsize(qh, qh->other_points); /* first centrum */
|
|
|
+ if (qh->DELAUNAY)
|
|
|
+ num--; /* drop last dimension */
|
|
|
+ qh_fprintf(qh, fp, 9061, "%d\n%d %d %d\n", num, qh->printoutvar
|
|
|
+ + numfacets - numsimplicial, numsimplicial + numridges, totneighbors/2);
|
|
|
+ }
|
|
|
+ FORALLpoints
|
|
|
+ qh_printpointid(qh, qh->fout, NULL, num, point, qh_IDunknown);
|
|
|
+ FOREACHpoint_(qh->other_points)
|
|
|
+ qh_printpointid(qh, qh->fout, NULL, num, point, qh_IDunknown);
|
|
|
+ if (format == qh_PRINTtriangles && qh->hull_dim > 2) {
|
|
|
+ FORALLfacets {
|
|
|
+ if (!facet->simplicial && facet->visitid)
|
|
|
+ qh_printcenter(qh, qh->fout, format, NULL, facet);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ break;
|
|
|
+ case qh_PRINTpointnearest:
|
|
|
+ qh_fprintf(qh, fp, 9062, "%d\n", numcoplanars);
|
|
|
+ break;
|
|
|
+ case qh_PRINTpoints:
|
|
|
+ if (!qh->VORONOI)
|
|
|
+ goto LABELnoformat;
|
|
|
+ if (qh->CDDoutput)
|
|
|
+ qh_fprintf(qh, fp, 9063, "%s | %s\nbegin\n%d %d real\n", qh->rbox_command,
|
|
|
+ qh->qhull_command, numfacets, qh->hull_dim);
|
|
|
+ else
|
|
|
+ qh_fprintf(qh, fp, 9064, "%d\n%d\n", qh->hull_dim-1, numfacets);
|
|
|
+ break;
|
|
|
+ case qh_PRINTvertices:
|
|
|
+ qh_fprintf(qh, fp, 9065, "%d\n", numfacets);
|
|
|
+ break;
|
|
|
+ case qh_PRINTsummary:
|
|
|
+ default:
|
|
|
+ LABELnoformat:
|
|
|
+ qh_fprintf(qh, qh->ferr, 6068, "qhull internal error (qh_printbegin): can not use this format for dimension %d\n",
|
|
|
+ qh->hull_dim);
|
|
|
+ qh_errexit(qh, qh_ERRqhull, NULL, NULL);
|
|
|
+ }
|
|
|
+} /* printbegin */
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printcenter">-</a>
|
|
|
+
|
|
|
+ qh_printcenter(qh, 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
|
|
|
+ Same as QhullFacet::printCenter
|
|
|
+*/
|
|
|
+void qh_printcenter(qhT *qh, FILE *fp, qh_PRINT format, const char *string, facetT *facet) {
|
|
|
+ int k, num;
|
|
|
+
|
|
|
+ if (qh->CENTERtype != qh_ASvoronoi && qh->CENTERtype != qh_AScentrum)
|
|
|
+ return;
|
|
|
+ if (string)
|
|
|
+ qh_fprintf(qh, fp, 9066, string);
|
|
|
+ if (qh->CENTERtype == qh_ASvoronoi) {
|
|
|
+ num= qh->hull_dim-1;
|
|
|
+ if (!facet->normal || !facet->upperdelaunay || !qh->ATinfinity) {
|
|
|
+ if (!facet->center)
|
|
|
+ facet->center= qh_facetcenter(qh, facet->vertices);
|
|
|
+ for (k=0; k < num; k++)
|
|
|
+ qh_fprintf(qh, fp, 9067, qh_REAL_1, facet->center[k]);
|
|
|
+ }else {
|
|
|
+ for (k=0; k < num; k++)
|
|
|
+ qh_fprintf(qh, fp, 9068, 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(qh, facet);
|
|
|
+ for (k=0; k < num; k++)
|
|
|
+ qh_fprintf(qh, fp, 9069, qh_REAL_1, facet->center[k]);
|
|
|
+ }
|
|
|
+ if (format == qh_PRINTgeom && num == 2)
|
|
|
+ qh_fprintf(qh, fp, 9070, " 0\n");
|
|
|
+ else
|
|
|
+ qh_fprintf(qh, fp, 9071, "\n");
|
|
|
+} /* printcenter */
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printcentrum">-</a>
|
|
|
+
|
|
|
+ qh_printcentrum(qh, 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(qhT *qh, 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(qh, facet);
|
|
|
+ centrum= facet->center;
|
|
|
+ }else {
|
|
|
+ centrum= qh_getcentrum(qh, facet);
|
|
|
+ tempcentrum= True;
|
|
|
+ }
|
|
|
+ qh_fprintf(qh, fp, 9072, "{appearance {-normal -edge normscale 0} ");
|
|
|
+ if (qh->firstcentrum) {
|
|
|
+ qh->firstcentrum= False;
|
|
|
+ qh_fprintf(qh, fp, 9073, "{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
|
|
|
+ qh_fprintf(qh, fp, 9074, "{INST geom { : centrum } transform { # f%d\n", facet->id);
|
|
|
+ apex= SETfirstt_(facet->vertices, vertexT);
|
|
|
+ qh_distplane(qh, apex->point, facet, &dist);
|
|
|
+ projpt= qh_projectpoint(qh, 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(qh, xaxis, xaxis);
|
|
|
+ qh_projectdim3(qh, normal, normal);
|
|
|
+ qh_normalize2(qh, normal, qh->PRINTdim, True, NULL, NULL);
|
|
|
+ }
|
|
|
+ qh_crossproduct(3, xaxis, normal, yaxis);
|
|
|
+ qh_fprintf(qh, fp, 9075, "%8.4g %8.4g %8.4g 0\n", xaxis[0], xaxis[1], xaxis[2]);
|
|
|
+ qh_fprintf(qh, fp, 9076, "%8.4g %8.4g %8.4g 0\n", yaxis[0], yaxis[1], yaxis[2]);
|
|
|
+ qh_fprintf(qh, fp, 9077, "%8.4g %8.4g %8.4g 0\n", normal[0], normal[1], normal[2]);
|
|
|
+ qh_printpoint3(qh, fp, centrum);
|
|
|
+ qh_fprintf(qh, fp, 9078, "1 }}}\n");
|
|
|
+ qh_memfree(qh, projpt, qh->normal_size);
|
|
|
+ qh_printpointvect(qh, fp, centrum, facet->normal, NULL, radius, green);
|
|
|
+ if (tempcentrum)
|
|
|
+ qh_memfree(qh, centrum, qh->normal_size);
|
|
|
+} /* printcentrum */
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printend">-</a>
|
|
|
+
|
|
|
+ qh_printend(qh, fp, format )
|
|
|
+ prints trailer for all output formats
|
|
|
+
|
|
|
+ see:
|
|
|
+ qh_printbegin() and qh_printafacet()
|
|
|
+
|
|
|
+*/
|
|
|
+void qh_printend(qhT *qh, FILE *fp, qh_PRINT format, facetT *facetlist, setT *facets, boolT printall) {
|
|
|
+ int num;
|
|
|
+ facetT *facet, **facetp;
|
|
|
+
|
|
|
+ if (!qh->printoutnum)
|
|
|
+ qh_fprintf(qh, qh->ferr, 7055, "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(qh, fp, facet,&num, printall);
|
|
|
+ FOREACHfacet_(facets)
|
|
|
+ qh_printend4geom(qh, fp, facet, &num, printall);
|
|
|
+ if (num != qh->ridgeoutnum || qh->printoutvar != qh->ridgeoutnum) {
|
|
|
+ qh_fprintf(qh, qh->ferr, 6069, "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, qh_ERRqhull, NULL, NULL);
|
|
|
+ }
|
|
|
+ }else
|
|
|
+ qh_fprintf(qh, fp, 9079, "}\n");
|
|
|
+ break;
|
|
|
+ case qh_PRINTinner:
|
|
|
+ case qh_PRINTnormals:
|
|
|
+ case qh_PRINTouter:
|
|
|
+ if (qh->CDDoutput)
|
|
|
+ qh_fprintf(qh, fp, 9080, "end\n");
|
|
|
+ break;
|
|
|
+ case qh_PRINTmaple:
|
|
|
+ qh_fprintf(qh, fp, 9081, "));\n");
|
|
|
+ break;
|
|
|
+ case qh_PRINTmathematica:
|
|
|
+ qh_fprintf(qh, fp, 9082, "}\n");
|
|
|
+ break;
|
|
|
+ case qh_PRINTpoints:
|
|
|
+ if (qh->CDDoutput)
|
|
|
+ qh_fprintf(qh, fp, 9083, "end\n");
|
|
|
+ break;
|
|
|
+ default:
|
|
|
+ break;
|
|
|
+ }
|
|
|
+} /* printend */
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printend4geom">-</a>
|
|
|
+
|
|
|
+ qh_printend4geom(qh, 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(qhT *qh, 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(qh, 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)
|
|
|
+ qh_fprintf(qh, fp, 9084, "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)
|
|
|
+ qh_fprintf(qh, fp, 9085, "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_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printextremes">-</a>
|
|
|
+
|
|
|
+ qh_printextremes(qh, 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(qhT *qh, FILE *fp, facetT *facetlist, setT *facets, boolT 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, qh->other_points);
|
|
|
+
|
|
|
+ points= qh_settemp(qh, allpoints);
|
|
|
+ qh_setzero(qh, points, 0, allpoints);
|
|
|
+ vertices= qh_facetvertices(qh, facetlist, facets, printall);
|
|
|
+ FOREACHvertex_(vertices) {
|
|
|
+ id= qh_pointid(qh, vertex->point);
|
|
|
+ if (id >= 0) {
|
|
|
+ SETelem_(points, id)= vertex->point;
|
|
|
+ numpoints++;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ qh_settempfree(qh, &vertices);
|
|
|
+ qh_fprintf(qh, fp, 9086, "%d\n", numpoints);
|
|
|
+ FOREACHpoint_i_(qh, points) {
|
|
|
+ if (point)
|
|
|
+ qh_fprintf(qh, fp, 9087, "%d\n", point_i);
|
|
|
+ }
|
|
|
+ qh_settempfree(qh, &points);
|
|
|
+} /* printextremes */
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printextremes_2d">-</a>
|
|
|
+
|
|
|
+ qh_printextremes_2d(qh, 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(qhT *qh, FILE *fp, facetT *facetlist, setT *facets, boolT printall) {
|
|
|
+ int numfacets, numridges, totneighbors, numcoplanars, numsimplicial, numtricoplanars;
|
|
|
+ setT *vertices;
|
|
|
+ facetT *facet, *startfacet, *nextfacet;
|
|
|
+ vertexT *vertexA, *vertexB;
|
|
|
+
|
|
|
+ qh_countfacets(qh, facetlist, facets, printall, &numfacets, &numsimplicial,
|
|
|
+ &totneighbors, &numridges, &numcoplanars, &numtricoplanars); /* marks qh->visit_id */
|
|
|
+ vertices= qh_facetvertices(qh, facetlist, facets, printall);
|
|
|
+ qh_fprintf(qh, fp, 9088, "%d\n", qh_setsize(qh, vertices));
|
|
|
+ qh_settempfree(qh, &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) {
|
|
|
+ qh_fprintf(qh, qh->ferr, 6218, "qhull internal error (qh_printextremes_2d): loop in facet list. facet %d nextfacet %d\n",
|
|
|
+ facet->id, nextfacet->id);
|
|
|
+ qh_errexit2(qh, qh_ERRqhull, facet, nextfacet);
|
|
|
+ }
|
|
|
+ if (facet->visitid) {
|
|
|
+ if (vertexA->visitid != qh->vertex_visit) {
|
|
|
+ vertexA->visitid= qh->vertex_visit;
|
|
|
+ qh_fprintf(qh, fp, 9089, "%d\n", qh_pointid(qh, vertexA->point));
|
|
|
+ }
|
|
|
+ if (vertexB->visitid != qh->vertex_visit) {
|
|
|
+ vertexB->visitid= qh->vertex_visit;
|
|
|
+ qh_fprintf(qh, fp, 9090, "%d\n", qh_pointid(qh, vertexB->point));
|
|
|
+ }
|
|
|
+ }
|
|
|
+ facet->visitid= qh->visit_id;
|
|
|
+ facet= nextfacet;
|
|
|
+ }while (facet && facet != startfacet);
|
|
|
+} /* printextremes_2d */
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printextremes_d">-</a>
|
|
|
+
|
|
|
+ qh_printextremes_d(qh, 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(qhT *qh, FILE *fp, facetT *facetlist, setT *facets, boolT printall) {
|
|
|
+ setT *vertices;
|
|
|
+ vertexT *vertex, **vertexp;
|
|
|
+ boolT upperseen, lowerseen;
|
|
|
+ facetT *neighbor, **neighborp;
|
|
|
+ int numpoints=0;
|
|
|
+
|
|
|
+ vertices= qh_facetvertices(qh, facetlist, facets, printall);
|
|
|
+ qh_vertexneighbors(qh);
|
|
|
+ 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;
|
|
|
+ }
|
|
|
+ qh_fprintf(qh, fp, 9091, "%d\n", numpoints);
|
|
|
+ FOREACHvertex_(vertices) {
|
|
|
+ if (vertex->seen)
|
|
|
+ qh_fprintf(qh, fp, 9092, "%d\n", qh_pointid(qh, vertex->point));
|
|
|
+ }
|
|
|
+ qh_settempfree(qh, &vertices);
|
|
|
+} /* printextremes_d */
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printfacet">-</a>
|
|
|
+
|
|
|
+ qh_printfacet(qh, fp, facet )
|
|
|
+ prints all fields of a facet to fp
|
|
|
+
|
|
|
+ notes:
|
|
|
+ ridges printed in neighbor order
|
|
|
+*/
|
|
|
+void qh_printfacet(qhT *qh, FILE *fp, facetT *facet) {
|
|
|
+
|
|
|
+ qh_printfacetheader(qh, fp, facet);
|
|
|
+ if (facet->ridges)
|
|
|
+ qh_printfacetridges(qh, fp, facet);
|
|
|
+} /* printfacet */
|
|
|
+
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printfacet2geom">-</a>
|
|
|
+
|
|
|
+ qh_printfacet2geom(qh, fp, facet, color )
|
|
|
+ print facet as part of a 2-d VECT for Geomview
|
|
|
+
|
|
|
+ notes:
|
|
|
+ assume precise calculations in io_r.c with roundoff covered by qh_GEOMepsilon
|
|
|
+ mindist is calculated within io_r.c. maxoutside is calculated elsewhere
|
|
|
+ so a DISTround error may have occurred.
|
|
|
+*/
|
|
|
+void qh_printfacet2geom(qhT *qh, FILE *fp, facetT *facet, realT color[3]) {
|
|
|
+ pointT *point0, *point1;
|
|
|
+ realT mindist, innerplane, outerplane;
|
|
|
+ int k;
|
|
|
+
|
|
|
+ qh_facet2point(qh, facet, &point0, &point1, &mindist);
|
|
|
+ qh_geomplanes(qh, facet, &outerplane, &innerplane);
|
|
|
+ if (qh->PRINTouter || (!qh->PRINTnoplanes && !qh->PRINTinner))
|
|
|
+ qh_printfacet2geom_points(qh, 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(qh, fp, point0, point1, facet, innerplane, color);
|
|
|
+ }
|
|
|
+ qh_memfree(qh, point1, qh->normal_size);
|
|
|
+ qh_memfree(qh, point0, qh->normal_size);
|
|
|
+} /* printfacet2geom */
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printfacet2geom_points">-</a>
|
|
|
+
|
|
|
+ qh_printfacet2geom_points(qh, 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(qhT *qh, FILE *fp, pointT *point1, pointT *point2,
|
|
|
+ facetT *facet, realT offset, realT color[3]) {
|
|
|
+ pointT *p1= point1, *p2= point2;
|
|
|
+
|
|
|
+ qh_fprintf(qh, fp, 9093, "VECT 1 2 1 2 1 # f%d\n", facet->id);
|
|
|
+ if (offset != 0.0) {
|
|
|
+ p1= qh_projectpoint(qh, p1, facet, -offset);
|
|
|
+ p2= qh_projectpoint(qh, p2, facet, -offset);
|
|
|
+ }
|
|
|
+ qh_fprintf(qh, fp, 9094, "%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(qh, p1, qh->normal_size);
|
|
|
+ qh_memfree(qh, p2, qh->normal_size);
|
|
|
+ }
|
|
|
+ qh_fprintf(qh, fp, 9095, "%8.4g %8.4g %8.4g 1.0\n", color[0], color[1], color[2]);
|
|
|
+} /* printfacet2geom_points */
|
|
|
+
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printfacet2math">-</a>
|
|
|
+
|
|
|
+ qh_printfacet2math(qh, fp, facet, format, notfirst )
|
|
|
+ print 2-d Maple or Mathematica output for a facet
|
|
|
+ may be non-simplicial
|
|
|
+
|
|
|
+ notes:
|
|
|
+ use %16.8f since Mathematica 2.2 does not handle exponential format
|
|
|
+ see qh_printfacet3math
|
|
|
+*/
|
|
|
+void qh_printfacet2math(qhT *qh, FILE *fp, facetT *facet, qh_PRINT format, int notfirst) {
|
|
|
+ pointT *point0, *point1;
|
|
|
+ realT mindist;
|
|
|
+ const char *pointfmt;
|
|
|
+
|
|
|
+ qh_facet2point(qh, facet, &point0, &point1, &mindist);
|
|
|
+ if (notfirst)
|
|
|
+ qh_fprintf(qh, fp, 9096, ",");
|
|
|
+ if (format == qh_PRINTmaple)
|
|
|
+ pointfmt= "[[%16.8f, %16.8f], [%16.8f, %16.8f]]\n";
|
|
|
+ else
|
|
|
+ pointfmt= "Line[{{%16.8f, %16.8f}, {%16.8f, %16.8f}}]\n";
|
|
|
+ qh_fprintf(qh, fp, 9097, pointfmt, point0[0], point0[1], point1[0], point1[1]);
|
|
|
+ qh_memfree(qh, point1, qh->normal_size);
|
|
|
+ qh_memfree(qh, point0, qh->normal_size);
|
|
|
+} /* printfacet2math */
|
|
|
+
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printfacet3geom_nonsimplicial">-</a>
|
|
|
+
|
|
|
+ qh_printfacet3geom_nonsimplicial(qh, 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(qhT *qh, 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(qh, facet, &outerplane, &innerplane);
|
|
|
+ vertices= qh_facet3vertex(qh, facet); /* oriented */
|
|
|
+ cntvertices= qh_setsize(qh, vertices);
|
|
|
+ projectedpoints= qh_settemp(qh, cntvertices);
|
|
|
+ FOREACHvertex_(vertices) {
|
|
|
+ zinc_(Zdistio);
|
|
|
+ qh_distplane(qh, vertex->point, facet, &dist);
|
|
|
+ projpt= qh_projectpoint(qh, vertex->point, facet, dist);
|
|
|
+ qh_setappend(qh, &projectedpoints, projpt);
|
|
|
+ }
|
|
|
+ if (qh->PRINTouter || (!qh->PRINTnoplanes && !qh->PRINTinner))
|
|
|
+ qh_printfacet3geom_points(qh, 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(qh, fp, projectedpoints, facet, innerplane, color);
|
|
|
+ }
|
|
|
+ FOREACHpoint_(projectedpoints)
|
|
|
+ qh_memfree(qh, point, qh->normal_size);
|
|
|
+ qh_settempfree(qh, &projectedpoints);
|
|
|
+ qh_settempfree(qh, &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(qh, fp, facet, neighbor, ridge->vertices, black);
|
|
|
+ if (qh->PRINTridges) {
|
|
|
+ vertexA= SETfirstt_(ridge->vertices, vertexT);
|
|
|
+ vertexB= SETsecondt_(ridge->vertices, vertexT);
|
|
|
+ qh_printline3geom(qh, fp, vertexA->point, vertexB->point, green);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+} /* printfacet3geom_nonsimplicial */
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printfacet3geom_points">-</a>
|
|
|
+
|
|
|
+ qh_printfacet3geom_points(qh, 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(qhT *qh, FILE *fp, setT *points, facetT *facet, realT offset, realT color[3]) {
|
|
|
+ int k, n= qh_setsize(qh, points), i;
|
|
|
+ pointT *point, **pointp;
|
|
|
+ setT *printpoints;
|
|
|
+
|
|
|
+ qh_fprintf(qh, fp, 9098, "{ OFF %d 1 1 # f%d\n", n, facet->id);
|
|
|
+ if (offset != 0.0) {
|
|
|
+ printpoints= qh_settemp(qh, n);
|
|
|
+ FOREACHpoint_(points)
|
|
|
+ qh_setappend(qh, &printpoints, qh_projectpoint(qh, point, facet, -offset));
|
|
|
+ }else
|
|
|
+ printpoints= points;
|
|
|
+ FOREACHpoint_(printpoints) {
|
|
|
+ for (k=0; k < qh->hull_dim; k++) {
|
|
|
+ if (k == qh->DROPdim)
|
|
|
+ qh_fprintf(qh, fp, 9099, "0 ");
|
|
|
+ else
|
|
|
+ qh_fprintf(qh, fp, 9100, "%8.4g ", point[k]);
|
|
|
+ }
|
|
|
+ if (printpoints != points)
|
|
|
+ qh_memfree(qh, point, qh->normal_size);
|
|
|
+ qh_fprintf(qh, fp, 9101, "\n");
|
|
|
+ }
|
|
|
+ if (printpoints != points)
|
|
|
+ qh_settempfree(qh, &printpoints);
|
|
|
+ qh_fprintf(qh, fp, 9102, "%d ", n);
|
|
|
+ for (i=0; i < n; i++)
|
|
|
+ qh_fprintf(qh, fp, 9103, "%d ", i);
|
|
|
+ qh_fprintf(qh, fp, 9104, "%8.4g %8.4g %8.4g 1.0 }\n", color[0], color[1], color[2]);
|
|
|
+} /* printfacet3geom_points */
|
|
|
+
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printfacet3geom_simplicial">-</a>
|
|
|
+
|
|
|
+ qh_printfacet3geom_simplicial(qh )
|
|
|
+ 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_r.c with roundoff covered by qh_GEOMepsilon
|
|
|
+ innerplane may be off by qh->DISTround. Maxoutside is calculated elsewhere
|
|
|
+ so a DISTround error may have occurred.
|
|
|
+*/
|
|
|
+void qh_printfacet3geom_simplicial(qhT *qh, 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(qh, facet, &outerplane, &innerplane);
|
|
|
+ vertices= qh_facet3vertex(qh, facet);
|
|
|
+ points= qh_settemp(qh, qh->TEMPsize);
|
|
|
+ FOREACHvertex_(vertices)
|
|
|
+ qh_setappend(qh, &points, vertex->point);
|
|
|
+ if (qh->PRINTouter || (!qh->PRINTnoplanes && !qh->PRINTinner))
|
|
|
+ qh_printfacet3geom_points(qh, 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(qh, fp, points, facet, innerplane, color);
|
|
|
+ }
|
|
|
+ qh_settempfree(qh, &points);
|
|
|
+ qh_settempfree(qh, &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(qh, facet->vertices, qh->hull_dim,
|
|
|
+ SETindex_(facet->neighbors, neighbor), 0);
|
|
|
+ if (qh->DOintersections)
|
|
|
+ qh_printhyperplaneintersection(qh, fp, facet, neighbor, vertices, black);
|
|
|
+ if (qh->PRINTridges) {
|
|
|
+ vertexA= SETfirstt_(vertices, vertexT);
|
|
|
+ vertexB= SETsecondt_(vertices, vertexT);
|
|
|
+ qh_printline3geom(qh, fp, vertexA->point, vertexB->point, green);
|
|
|
+ }
|
|
|
+ qh_setfree(qh, &vertices);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+} /* printfacet3geom_simplicial */
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printfacet3math">-</a>
|
|
|
+
|
|
|
+ qh_printfacet3math(qh, fp, facet, notfirst )
|
|
|
+ print 3-d Maple or Mathematica output for a facet
|
|
|
+
|
|
|
+ notes:
|
|
|
+ may be non-simplicial
|
|
|
+ use %16.8f since Mathematica 2.2 does not handle exponential format
|
|
|
+ see qh_printfacet2math
|
|
|
+*/
|
|
|
+void qh_printfacet3math(qhT *qh, FILE *fp, facetT *facet, qh_PRINT format, int notfirst) {
|
|
|
+ vertexT *vertex, **vertexp;
|
|
|
+ setT *points, *vertices;
|
|
|
+ pointT *point, **pointp;
|
|
|
+ boolT firstpoint= True;
|
|
|
+ realT dist;
|
|
|
+ const char *pointfmt, *endfmt;
|
|
|
+
|
|
|
+ if (notfirst)
|
|
|
+ qh_fprintf(qh, fp, 9105, ",\n");
|
|
|
+ vertices= qh_facet3vertex(qh, facet);
|
|
|
+ points= qh_settemp(qh, qh_setsize(qh, vertices));
|
|
|
+ FOREACHvertex_(vertices) {
|
|
|
+ zinc_(Zdistio);
|
|
|
+ qh_distplane(qh, vertex->point, facet, &dist);
|
|
|
+ point= qh_projectpoint(qh, vertex->point, facet, dist);
|
|
|
+ qh_setappend(qh, &points, point);
|
|
|
+ }
|
|
|
+ if (format == qh_PRINTmaple) {
|
|
|
+ qh_fprintf(qh, fp, 9106, "[");
|
|
|
+ pointfmt= "[%16.8f, %16.8f, %16.8f]";
|
|
|
+ endfmt= "]";
|
|
|
+ }else {
|
|
|
+ qh_fprintf(qh, fp, 9107, "Polygon[{");
|
|
|
+ pointfmt= "{%16.8f, %16.8f, %16.8f}";
|
|
|
+ endfmt= "}]";
|
|
|
+ }
|
|
|
+ FOREACHpoint_(points) {
|
|
|
+ if (firstpoint)
|
|
|
+ firstpoint= False;
|
|
|
+ else
|
|
|
+ qh_fprintf(qh, fp, 9108, ",\n");
|
|
|
+ qh_fprintf(qh, fp, 9109, pointfmt, point[0], point[1], point[2]);
|
|
|
+ }
|
|
|
+ FOREACHpoint_(points)
|
|
|
+ qh_memfree(qh, point, qh->normal_size);
|
|
|
+ qh_settempfree(qh, &points);
|
|
|
+ qh_settempfree(qh, &vertices);
|
|
|
+ qh_fprintf(qh, fp, 9110, "%s", endfmt);
|
|
|
+} /* printfacet3math */
|
|
|
+
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printfacet3vertex">-</a>
|
|
|
+
|
|
|
+ qh_printfacet3vertex(qh, 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(qhT *qh, FILE *fp, facetT *facet, qh_PRINT format) {
|
|
|
+ vertexT *vertex, **vertexp;
|
|
|
+ setT *vertices;
|
|
|
+
|
|
|
+ vertices= qh_facet3vertex(qh, facet);
|
|
|
+ if (format == qh_PRINToff)
|
|
|
+ qh_fprintf(qh, fp, 9111, "%d ", qh_setsize(qh, vertices));
|
|
|
+ FOREACHvertex_(vertices)
|
|
|
+ qh_fprintf(qh, fp, 9112, "%d ", qh_pointid(qh, vertex->point));
|
|
|
+ qh_fprintf(qh, fp, 9113, "\n");
|
|
|
+ qh_settempfree(qh, &vertices);
|
|
|
+} /* printfacet3vertex */
|
|
|
+
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printfacet4geom_nonsimplicial">-</a>
|
|
|
+
|
|
|
+ qh_printfacet4geom_nonsimplicial(qh )
|
|
|
+ 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(qhT *qh, 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(qh, fp, facet, neighbor, ridge->vertices, color);
|
|
|
+ else {
|
|
|
+ if (qh->DROPdim >= 0)
|
|
|
+ qh_fprintf(qh, fp, 9114, "OFF 3 1 1 # f%d\n", facet->id);
|
|
|
+ else {
|
|
|
+ qh->printoutvar++;
|
|
|
+ qh_fprintf(qh, fp, 9115, "# r%d between f%d f%d\n", ridge->id, facet->id, neighbor->id);
|
|
|
+ }
|
|
|
+ FOREACHvertex_(ridge->vertices) {
|
|
|
+ zinc_(Zdistio);
|
|
|
+ qh_distplane(qh, vertex->point,facet, &dist);
|
|
|
+ point=qh_projectpoint(qh, vertex->point,facet, dist);
|
|
|
+ for (k=0; k < qh->hull_dim; k++) {
|
|
|
+ if (k != qh->DROPdim)
|
|
|
+ qh_fprintf(qh, fp, 9116, "%8.4g ", point[k]);
|
|
|
+ }
|
|
|
+ qh_fprintf(qh, fp, 9117, "\n");
|
|
|
+ qh_memfree(qh, point, qh->normal_size);
|
|
|
+ }
|
|
|
+ if (qh->DROPdim >= 0)
|
|
|
+ qh_fprintf(qh, fp, 9118, "3 0 1 2 %8.4g %8.4g %8.4g\n", color[0], color[1], color[2]);
|
|
|
+ }
|
|
|
+ }
|
|
|
+} /* printfacet4geom_nonsimplicial */
|
|
|
+
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printfacet4geom_simplicial">-</a>
|
|
|
+
|
|
|
+ qh_printfacet4geom_simplicial(qh, 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(qhT *qh, 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(qh, facet->vertices, qh->hull_dim,
|
|
|
+ SETindex_(facet->neighbors, neighbor), 0);
|
|
|
+ if (qh->DOintersections)
|
|
|
+ qh_printhyperplaneintersection(qh, fp, facet, neighbor, vertices, color);
|
|
|
+ else {
|
|
|
+ if (qh->DROPdim >= 0)
|
|
|
+ qh_fprintf(qh, fp, 9119, "OFF 3 1 1 # ridge between f%d f%d\n",
|
|
|
+ facet->id, neighbor->id);
|
|
|
+ else {
|
|
|
+ qh->printoutvar++;
|
|
|
+ qh_fprintf(qh, fp, 9120, "# 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)
|
|
|
+ qh_fprintf(qh, fp, 9121, "%8.4g ", vertex->point[k]);
|
|
|
+ }
|
|
|
+ qh_fprintf(qh, fp, 9122, "\n");
|
|
|
+ }
|
|
|
+ if (qh->DROPdim >= 0)
|
|
|
+ qh_fprintf(qh, fp, 9123, "3 0 1 2 %8.4g %8.4g %8.4g\n", color[0], color[1], color[2]);
|
|
|
+ }
|
|
|
+ qh_setfree(qh, &vertices);
|
|
|
+ }
|
|
|
+} /* printfacet4geom_simplicial */
|
|
|
+
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printfacetNvertex_nonsimplicial">-</a>
|
|
|
+
|
|
|
+ qh_printfacetNvertex_nonsimplicial(qh, fp, facet, id, format )
|
|
|
+ print vertices for an N-d non-simplicial facet
|
|
|
+ triangulates each ridge to the id
|
|
|
+*/
|
|
|
+void qh_printfacetNvertex_nonsimplicial(qhT *qh, FILE *fp, facetT *facet, int id, qh_PRINT format) {
|
|
|
+ vertexT *vertex, **vertexp;
|
|
|
+ ridgeT *ridge, **ridgep;
|
|
|
+
|
|
|
+ if (facet->visible && qh->NEWfacets)
|
|
|
+ return;
|
|
|
+ FOREACHridge_(facet->ridges) {
|
|
|
+ if (format == qh_PRINTtriangles)
|
|
|
+ qh_fprintf(qh, fp, 9124, "%d ", qh->hull_dim);
|
|
|
+ qh_fprintf(qh, fp, 9125, "%d ", id);
|
|
|
+ if ((ridge->top == facet) ^ qh_ORIENTclock) {
|
|
|
+ FOREACHvertex_(ridge->vertices)
|
|
|
+ qh_fprintf(qh, fp, 9126, "%d ", qh_pointid(qh, vertex->point));
|
|
|
+ }else {
|
|
|
+ FOREACHvertexreverse12_(ridge->vertices)
|
|
|
+ qh_fprintf(qh, fp, 9127, "%d ", qh_pointid(qh, vertex->point));
|
|
|
+ }
|
|
|
+ qh_fprintf(qh, fp, 9128, "\n");
|
|
|
+ }
|
|
|
+} /* printfacetNvertex_nonsimplicial */
|
|
|
+
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printfacetNvertex_simplicial">-</a>
|
|
|
+
|
|
|
+ qh_printfacetNvertex_simplicial(qh, 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(qhT *qh, FILE *fp, facetT *facet, qh_PRINT format) {
|
|
|
+ vertexT *vertex, **vertexp;
|
|
|
+
|
|
|
+ if (format == qh_PRINToff || format == qh_PRINTtriangles)
|
|
|
+ qh_fprintf(qh, fp, 9129, "%d ", qh_setsize(qh, facet->vertices));
|
|
|
+ if ((facet->toporient ^ qh_ORIENTclock)
|
|
|
+ || (qh->hull_dim > 2 && !facet->simplicial)) {
|
|
|
+ FOREACHvertex_(facet->vertices)
|
|
|
+ qh_fprintf(qh, fp, 9130, "%d ", qh_pointid(qh, vertex->point));
|
|
|
+ }else {
|
|
|
+ FOREACHvertexreverse12_(facet->vertices)
|
|
|
+ qh_fprintf(qh, fp, 9131, "%d ", qh_pointid(qh, vertex->point));
|
|
|
+ }
|
|
|
+ qh_fprintf(qh, fp, 9132, "\n");
|
|
|
+} /* printfacetNvertex_simplicial */
|
|
|
+
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printfacetheader">-</a>
|
|
|
+
|
|
|
+ qh_printfacetheader(qh, fp, facet )
|
|
|
+ prints header fields of a facet to fp
|
|
|
+
|
|
|
+ notes:
|
|
|
+ for 'f' output and debugging
|
|
|
+ Same as QhullFacet::printHeader()
|
|
|
+*/
|
|
|
+void qh_printfacetheader(qhT *qh, FILE *fp, facetT *facet) {
|
|
|
+ pointT *point, **pointp, *furthest;
|
|
|
+ facetT *neighbor, **neighborp;
|
|
|
+ realT dist;
|
|
|
+
|
|
|
+ if (facet == qh_MERGEridge) {
|
|
|
+ qh_fprintf(qh, fp, 9133, " MERGEridge\n");
|
|
|
+ return;
|
|
|
+ }else if (facet == qh_DUPLICATEridge) {
|
|
|
+ qh_fprintf(qh, fp, 9134, " DUPLICATEridge\n");
|
|
|
+ return;
|
|
|
+ }else if (!facet) {
|
|
|
+ qh_fprintf(qh, fp, 9135, " NULLfacet\n");
|
|
|
+ return;
|
|
|
+ }
|
|
|
+ qh->old_randomdist= qh->RANDOMdist;
|
|
|
+ qh->RANDOMdist= False;
|
|
|
+ qh_fprintf(qh, fp, 9136, "- f%d\n", facet->id);
|
|
|
+ qh_fprintf(qh, fp, 9137, " - flags:");
|
|
|
+ if (facet->toporient)
|
|
|
+ qh_fprintf(qh, fp, 9138, " top");
|
|
|
+ else
|
|
|
+ qh_fprintf(qh, fp, 9139, " bottom");
|
|
|
+ if (facet->simplicial)
|
|
|
+ qh_fprintf(qh, fp, 9140, " simplicial");
|
|
|
+ if (facet->tricoplanar)
|
|
|
+ qh_fprintf(qh, fp, 9141, " tricoplanar");
|
|
|
+ if (facet->upperdelaunay)
|
|
|
+ qh_fprintf(qh, fp, 9142, " upperDelaunay");
|
|
|
+ if (facet->visible)
|
|
|
+ qh_fprintf(qh, fp, 9143, " visible");
|
|
|
+ if (facet->newfacet)
|
|
|
+ qh_fprintf(qh, fp, 9144, " newfacet");
|
|
|
+ if (facet->tested)
|
|
|
+ qh_fprintf(qh, fp, 9145, " tested");
|
|
|
+ if (!facet->good)
|
|
|
+ qh_fprintf(qh, fp, 9146, " notG");
|
|
|
+ if (facet->seen && qh->IStracing)
|
|
|
+ qh_fprintf(qh, fp, 9147, " seen");
|
|
|
+ if (facet->seen2 && qh->IStracing)
|
|
|
+ qh_fprintf(qh, fp, 9418, " seen2");
|
|
|
+ if (facet->isarea)
|
|
|
+ qh_fprintf(qh, fp, 9419, " isarea");
|
|
|
+ if (facet->coplanarhorizon)
|
|
|
+ qh_fprintf(qh, fp, 9148, " coplanarhorizon");
|
|
|
+ if (facet->mergehorizon)
|
|
|
+ qh_fprintf(qh, fp, 9149, " mergehorizon");
|
|
|
+ if (facet->cycledone)
|
|
|
+ qh_fprintf(qh, fp, 9420, " cycledone");
|
|
|
+ if (facet->keepcentrum)
|
|
|
+ qh_fprintf(qh, fp, 9150, " keepcentrum");
|
|
|
+ if (facet->dupridge)
|
|
|
+ qh_fprintf(qh, fp, 9151, " dupridge");
|
|
|
+ if (facet->mergeridge && !facet->mergeridge2)
|
|
|
+ qh_fprintf(qh, fp, 9152, " mergeridge1");
|
|
|
+ if (facet->mergeridge2)
|
|
|
+ qh_fprintf(qh, fp, 9153, " mergeridge2");
|
|
|
+ if (facet->newmerge)
|
|
|
+ qh_fprintf(qh, fp, 9154, " newmerge");
|
|
|
+ if (facet->flipped)
|
|
|
+ qh_fprintf(qh, fp, 9155, " flipped");
|
|
|
+ if (facet->notfurthest)
|
|
|
+ qh_fprintf(qh, fp, 9156, " notfurthest");
|
|
|
+ if (facet->degenerate)
|
|
|
+ qh_fprintf(qh, fp, 9157, " degenerate");
|
|
|
+ if (facet->redundant)
|
|
|
+ qh_fprintf(qh, fp, 9158, " redundant");
|
|
|
+ qh_fprintf(qh, fp, 9159, "\n");
|
|
|
+ if (facet->isarea)
|
|
|
+ qh_fprintf(qh, fp, 9160, " - area: %2.2g\n", facet->f.area);
|
|
|
+ else if (qh->NEWfacets && facet->visible && facet->f.replace)
|
|
|
+ qh_fprintf(qh, fp, 9161, " - replacement: f%d\n", facet->f.replace->id);
|
|
|
+ else if (facet->newfacet) {
|
|
|
+ if (facet->f.samecycle && facet->f.samecycle != facet)
|
|
|
+ qh_fprintf(qh, fp, 9162, " - shares same visible/horizon as f%d\n", facet->f.samecycle->id);
|
|
|
+ }else if (facet->tricoplanar /* !isarea */) {
|
|
|
+ if (facet->f.triowner)
|
|
|
+ qh_fprintf(qh, fp, 9163, " - owner of normal & centrum is facet f%d\n", facet->f.triowner->id);
|
|
|
+ }else if (facet->f.newcycle)
|
|
|
+ qh_fprintf(qh, fp, 9164, " - was horizon to f%d\n", facet->f.newcycle->id);
|
|
|
+ if (facet->nummerge == qh_MAXnummerge)
|
|
|
+ qh_fprintf(qh, fp, 9427, " - merges: %dmax\n", qh_MAXnummerge);
|
|
|
+ else if (facet->nummerge)
|
|
|
+ qh_fprintf(qh, fp, 9165, " - merges: %d\n", facet->nummerge);
|
|
|
+ qh_printpointid(qh, fp, " - normal: ", qh->hull_dim, facet->normal, qh_IDunknown);
|
|
|
+ qh_fprintf(qh, fp, 9166, " - offset: %10.7g\n", facet->offset);
|
|
|
+ if (qh->CENTERtype == qh_ASvoronoi || facet->center)
|
|
|
+ qh_printcenter(qh, fp, qh_PRINTfacets, " - center: ", facet);
|
|
|
+#if qh_MAXoutside
|
|
|
+ if (facet->maxoutside > qh->DISTround) /* initial value */
|
|
|
+ qh_fprintf(qh, fp, 9167, " - maxoutside: %10.7g\n", facet->maxoutside);
|
|
|
+#endif
|
|
|
+ if (!SETempty_(facet->outsideset)) {
|
|
|
+ furthest= (pointT *)qh_setlast(facet->outsideset);
|
|
|
+ if (qh_setsize(qh, facet->outsideset) < 6) {
|
|
|
+ qh_fprintf(qh, fp, 9168, " - outside set(furthest p%d):\n", qh_pointid(qh, furthest));
|
|
|
+ FOREACHpoint_(facet->outsideset)
|
|
|
+ qh_printpoint(qh, fp, " ", point);
|
|
|
+ }else if (qh_setsize(qh, facet->outsideset) < 21) {
|
|
|
+ qh_printpoints(qh, fp, " - outside set:", facet->outsideset);
|
|
|
+ }else {
|
|
|
+ qh_fprintf(qh, fp, 9169, " - outside set: %d points.", qh_setsize(qh, facet->outsideset));
|
|
|
+ qh_printpoint(qh, fp, " Furthest", furthest);
|
|
|
+ }
|
|
|
+#if !qh_COMPUTEfurthest
|
|
|
+ qh_fprintf(qh, fp, 9170, " - furthest distance= %2.2g\n", facet->furthestdist);
|
|
|
+#endif
|
|
|
+ }
|
|
|
+ if (!SETempty_(facet->coplanarset)) {
|
|
|
+ furthest= (pointT *)qh_setlast(facet->coplanarset);
|
|
|
+ if (qh_setsize(qh, facet->coplanarset) < 6) {
|
|
|
+ qh_fprintf(qh, fp, 9171, " - coplanar set(furthest p%d):\n", qh_pointid(qh, furthest));
|
|
|
+ FOREACHpoint_(facet->coplanarset)
|
|
|
+ qh_printpoint(qh, fp, " ", point);
|
|
|
+ }else if (qh_setsize(qh, facet->coplanarset) < 21) {
|
|
|
+ qh_printpoints(qh, fp, " - coplanar set:", facet->coplanarset);
|
|
|
+ }else {
|
|
|
+ qh_fprintf(qh, fp, 9172, " - coplanar set: %d points.", qh_setsize(qh, facet->coplanarset));
|
|
|
+ qh_printpoint(qh, fp, " Furthest", furthest);
|
|
|
+ }
|
|
|
+ zinc_(Zdistio);
|
|
|
+ qh_distplane(qh, furthest, facet, &dist);
|
|
|
+ qh_fprintf(qh, fp, 9173, " furthest distance= %2.2g\n", dist);
|
|
|
+ }
|
|
|
+ qh_printvertices(qh, fp, " - vertices:", facet->vertices);
|
|
|
+ qh_fprintf(qh, fp, 9174, " - neighboring facets:");
|
|
|
+ FOREACHneighbor_(facet) {
|
|
|
+ if (neighbor == qh_MERGEridge)
|
|
|
+ qh_fprintf(qh, fp, 9175, " MERGEridge");
|
|
|
+ else if (neighbor == qh_DUPLICATEridge)
|
|
|
+ qh_fprintf(qh, fp, 9176, " DUPLICATEridge");
|
|
|
+ else
|
|
|
+ qh_fprintf(qh, fp, 9177, " f%d", neighbor->id);
|
|
|
+ }
|
|
|
+ qh_fprintf(qh, fp, 9178, "\n");
|
|
|
+ qh->RANDOMdist= qh->old_randomdist;
|
|
|
+} /* printfacetheader */
|
|
|
+
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printfacetridges">-</a>
|
|
|
+
|
|
|
+ qh_printfacetridges(qh, fp, facet )
|
|
|
+ prints ridges of a facet to fp
|
|
|
+
|
|
|
+ notes:
|
|
|
+ ridges printed in neighbor order
|
|
|
+ assumes the ridges exist
|
|
|
+ for 'f' output
|
|
|
+ same as QhullFacet::printRidges
|
|
|
+*/
|
|
|
+void qh_printfacetridges(qhT *qh, FILE *fp, facetT *facet) {
|
|
|
+ facetT *neighbor, **neighborp;
|
|
|
+ ridgeT *ridge, **ridgep;
|
|
|
+ int numridges= 0;
|
|
|
+ int n;
|
|
|
+
|
|
|
+ if (facet->visible && qh->NEWfacets) {
|
|
|
+ qh_fprintf(qh, fp, 9179, " - ridges (tentative ids):");
|
|
|
+ FOREACHridge_(facet->ridges)
|
|
|
+ qh_fprintf(qh, fp, 9180, " r%d", ridge->id);
|
|
|
+ qh_fprintf(qh, fp, 9181, "\n");
|
|
|
+ }else {
|
|
|
+ qh_fprintf(qh, fp, 9182, " - 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(qh, fp, ridge);
|
|
|
+ numridges++;
|
|
|
+ ridge= qh_nextridge3d(ridge, facet, NULL);
|
|
|
+ }
|
|
|
+ }else {
|
|
|
+ FOREACHneighbor_(facet) {
|
|
|
+ FOREACHridge_(facet->ridges) {
|
|
|
+ if (otherfacet_(ridge, facet) == neighbor && !ridge->seen) {
|
|
|
+ ridge->seen= True;
|
|
|
+ qh_printridge(qh, fp, ridge);
|
|
|
+ numridges++;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ n= qh_setsize(qh, facet->ridges);
|
|
|
+ if (n == 1 && facet->newfacet && qh->NEWtentative) {
|
|
|
+ qh_fprintf(qh, fp, 9411, " - horizon ridge to visible facet\n");
|
|
|
+ }
|
|
|
+ if (numridges != n) {
|
|
|
+ qh_fprintf(qh, fp, 9183, " - all ridges:");
|
|
|
+ FOREACHridge_(facet->ridges)
|
|
|
+ qh_fprintf(qh, fp, 9184, " r%d", ridge->id);
|
|
|
+ qh_fprintf(qh, fp, 9185, "\n");
|
|
|
+ }
|
|
|
+ /* non-3d ridges w/o non-simplicial neighbors */
|
|
|
+ FOREACHridge_(facet->ridges) {
|
|
|
+ if (!ridge->seen)
|
|
|
+ qh_printridge(qh, fp, ridge);
|
|
|
+ }
|
|
|
+ }
|
|
|
+} /* printfacetridges */
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printfacets">-</a>
|
|
|
+
|
|
|
+ qh_printfacets(qh, 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(qhT *qh, FILE *fp, qh_PRINT 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))
|
|
|
+ qh_fprintf(qh, qh->ferr, 7056, "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(qh, facetlist, facets, printall);
|
|
|
+ center= qh_getcenter(qh, vertices);
|
|
|
+ qh_fprintf(qh, fp, 9186, "%d 1\n", qh->hull_dim);
|
|
|
+ qh_printpointid(qh, fp, NULL, qh->hull_dim, center, qh_IDunknown);
|
|
|
+ qh_memfree(qh, center, qh->normal_size);
|
|
|
+ qh_settempfree(qh, &vertices);
|
|
|
+ }else if (format == qh_PRINTextremes) {
|
|
|
+ if (qh->DELAUNAY)
|
|
|
+ qh_printextremes_d(qh, fp, facetlist, facets, printall);
|
|
|
+ else if (qh->hull_dim == 2)
|
|
|
+ qh_printextremes_2d(qh, fp, facetlist, facets, printall);
|
|
|
+ else
|
|
|
+ qh_printextremes(qh, fp, facetlist, facets, printall);
|
|
|
+ }else if (format == qh_PRINToptions)
|
|
|
+ qh_fprintf(qh, fp, 9187, "Options selected for Qhull %s:\n%s\n", qh_version, qh->qhull_options);
|
|
|
+ else if (format == qh_PRINTpoints && !qh->VORONOI)
|
|
|
+ qh_printpoints_out(qh, fp, facetlist, facets, printall);
|
|
|
+ else if (format == qh_PRINTqhull)
|
|
|
+ qh_fprintf(qh, fp, 9188, "%s | %s\n", qh->rbox_command, qh->qhull_command);
|
|
|
+ else if (format == qh_PRINTsize) {
|
|
|
+ qh_fprintf(qh, fp, 9189, "0\n2 ");
|
|
|
+ qh_fprintf(qh, fp, 9190, qh_REAL_1, qh->totarea);
|
|
|
+ qh_fprintf(qh, fp, 9191, qh_REAL_1, qh->totvol);
|
|
|
+ qh_fprintf(qh, fp, 9192, "\n");
|
|
|
+ }else if (format == qh_PRINTsummary) {
|
|
|
+ qh_countfacets(qh, facetlist, facets, printall, &numfacets, &numsimplicial,
|
|
|
+ &totneighbors, &numridges, &numcoplanars, &numtricoplanars);
|
|
|
+ vertices= qh_facetvertices(qh, facetlist, facets, printall);
|
|
|
+ qh_fprintf(qh, fp, 9193, "10 %d %d %d %d %d %d %d %d %d %d\n2 ", qh->hull_dim,
|
|
|
+ qh->num_points + qh_setsize(qh, qh->other_points),
|
|
|
+ qh->num_vertices, qh->num_facets - qh->num_visible,
|
|
|
+ qh_setsize(qh, vertices), numfacets, numcoplanars,
|
|
|
+ numfacets - numsimplicial, zzval_(Zdelvertextot),
|
|
|
+ numtricoplanars);
|
|
|
+ qh_settempfree(qh, &vertices);
|
|
|
+ qh_outerinner(qh, NULL, &outerplane, &innerplane);
|
|
|
+ qh_fprintf(qh, fp, 9194, qh_REAL_2n, outerplane, innerplane);
|
|
|
+ }else if (format == qh_PRINTvneighbors)
|
|
|
+ qh_printvneighbors(qh, fp, facetlist, facets, printall);
|
|
|
+ else if (qh->VORONOI && format == qh_PRINToff)
|
|
|
+ qh_printvoronoi(qh, fp, format, facetlist, facets, printall);
|
|
|
+ else if (qh->VORONOI && format == qh_PRINTgeom) {
|
|
|
+ qh_printbegin(qh, fp, format, facetlist, facets, printall);
|
|
|
+ qh_printvoronoi(qh, fp, format, facetlist, facets, printall);
|
|
|
+ qh_printend(qh, fp, format, facetlist, facets, printall);
|
|
|
+ }else if (qh->VORONOI
|
|
|
+ && (format == qh_PRINTvertices || format == qh_PRINTinner || format == qh_PRINTouter))
|
|
|
+ qh_printvdiagram(qh, fp, format, facetlist, facets, printall);
|
|
|
+ else {
|
|
|
+ qh_printbegin(qh, fp, format, facetlist, facets, printall);
|
|
|
+ FORALLfacet_(facetlist)
|
|
|
+ qh_printafacet(qh, fp, format, facet, printall);
|
|
|
+ FOREACHfacet_(facets)
|
|
|
+ qh_printafacet(qh, fp, format, facet, printall);
|
|
|
+ qh_printend(qh, fp, format, facetlist, facets, printall);
|
|
|
+ }
|
|
|
+ qh->RANDOMdist= qh->old_randomdist;
|
|
|
+} /* printfacets */
|
|
|
+
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printhyperplaneintersection">-</a>
|
|
|
+
|
|
|
+ qh_printhyperplaneintersection(qh, 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(qhT *qh, 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(qh, facet1->normal, facet2->normal);
|
|
|
+ denominator= 1 - costheta * costheta;
|
|
|
+ i= qh_setsize(qh, vertices);
|
|
|
+ if (qh->hull_dim == 3)
|
|
|
+ qh_fprintf(qh, fp, 9195, "VECT 1 %d 1 %d 1 ", i, i);
|
|
|
+ else if (qh->hull_dim == 4 && qh->DROPdim >= 0)
|
|
|
+ qh_fprintf(qh, fp, 9196, "OFF 3 1 1 ");
|
|
|
+ else
|
|
|
+ qh->printoutvar++;
|
|
|
+ qh_fprintf(qh, fp, 9197, "# intersect f%d f%d\n", facet1->id, facet2->id);
|
|
|
+ mindenom= 1 / (10.0 * qh->MAXabs_coord);
|
|
|
+ FOREACHvertex_(vertices) {
|
|
|
+ zadd_(Zdistio, 2);
|
|
|
+ qh_distplane(qh, vertex->point, facet1, &dist1);
|
|
|
+ qh_distplane(qh, 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(qh, p, p);
|
|
|
+ qh_fprintf(qh, fp, 9198, "%8.4g %8.4g %8.4g # ", p[0], p[1], p[2]);
|
|
|
+ }else
|
|
|
+ qh_fprintf(qh, fp, 9199, "%8.4g %8.4g %8.4g %8.4g # ", p[0], p[1], p[2], p[3]);
|
|
|
+ if (nearzero1+nearzero2)
|
|
|
+ qh_fprintf(qh, fp, 9200, "p%d(coplanar facets)\n", qh_pointid(qh, vertex->point));
|
|
|
+ else
|
|
|
+ qh_fprintf(qh, fp, 9201, "projected p%d\n", qh_pointid(qh, vertex->point));
|
|
|
+ }
|
|
|
+ if (qh->hull_dim == 3)
|
|
|
+ qh_fprintf(qh, fp, 9202, "%8.4g %8.4g %8.4g 1.0\n", color[0], color[1], color[2]);
|
|
|
+ else if (qh->hull_dim == 4 && qh->DROPdim >= 0)
|
|
|
+ qh_fprintf(qh, fp, 9203, "3 0 1 2 %8.4g %8.4g %8.4g 1.0\n", color[0], color[1], color[2]);
|
|
|
+} /* printhyperplaneintersection */
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printline3geom">-</a>
|
|
|
+
|
|
|
+ qh_printline3geom(qh, 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(qhT *qh, FILE *fp, pointT *pointA, pointT *pointB, realT color[3]) {
|
|
|
+ int k;
|
|
|
+ realT pA[4], pB[4];
|
|
|
+
|
|
|
+ qh_projectdim3(qh, pointA, pA);
|
|
|
+ qh_projectdim3(qh, pointB, pB);
|
|
|
+ if ((fabs(pA[0] - pB[0]) > 1e-3) ||
|
|
|
+ (fabs(pA[1] - pB[1]) > 1e-3) ||
|
|
|
+ (fabs(pA[2] - pB[2]) > 1e-3)) {
|
|
|
+ qh_fprintf(qh, fp, 9204, "VECT 1 2 1 2 1\n");
|
|
|
+ for (k=0; k < 3; k++)
|
|
|
+ qh_fprintf(qh, fp, 9205, "%8.4g ", pB[k]);
|
|
|
+ qh_fprintf(qh, fp, 9206, " # p%d\n", qh_pointid(qh, pointB));
|
|
|
+ }else
|
|
|
+ qh_fprintf(qh, fp, 9207, "VECT 1 1 1 1 1\n");
|
|
|
+ for (k=0; k < 3; k++)
|
|
|
+ qh_fprintf(qh, fp, 9208, "%8.4g ", pA[k]);
|
|
|
+ qh_fprintf(qh, fp, 9209, " # p%d\n", qh_pointid(qh, pointA));
|
|
|
+ qh_fprintf(qh, fp, 9210, "%8.4g %8.4g %8.4g 1\n", color[0], color[1], color[2]);
|
|
|
+}
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printneighborhood">-</a>
|
|
|
+
|
|
|
+ qh_printneighborhood(qh, fp, format, facetA, facetB, printall )
|
|
|
+ print neighborhood of one or two facets
|
|
|
+
|
|
|
+ notes:
|
|
|
+ calls qh_findgood_all()
|
|
|
+ bumps qh.visit_id
|
|
|
+*/
|
|
|
+void qh_printneighborhood(qhT *qh, FILE *fp, qh_PRINT format, facetT *facetA, facetT *facetB, boolT printall) {
|
|
|
+ facetT *neighbor, **neighborp, *facet;
|
|
|
+ setT *facets;
|
|
|
+
|
|
|
+ if (format == qh_PRINTnone)
|
|
|
+ return;
|
|
|
+ qh_findgood_all(qh, qh->facet_list);
|
|
|
+ if (facetA == facetB)
|
|
|
+ facetB= NULL;
|
|
|
+ facets= qh_settemp(qh, 2*(qh_setsize(qh, 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(qh, &facets, facet);
|
|
|
+ }
|
|
|
+ FOREACHneighbor_(facet) {
|
|
|
+ if (neighbor->visitid == qh->visit_id)
|
|
|
+ continue;
|
|
|
+ neighbor->visitid= qh->visit_id;
|
|
|
+ if (printall || !qh_skipfacet(qh, neighbor))
|
|
|
+ qh_setappend(qh, &facets, neighbor);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ qh_printfacets(qh, fp, format, NULL, facets, printall);
|
|
|
+ qh_settempfree(qh, &facets);
|
|
|
+} /* printneighborhood */
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printpoint">-</a>
|
|
|
+
|
|
|
+ qh_printpoint(qh, fp, string, point )
|
|
|
+ qh_printpointid(qh, 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=qh_IDunknown(-1) or qh_IDnone(-3)
|
|
|
+
|
|
|
+ notes:
|
|
|
+ nop if point is NULL
|
|
|
+ Same as QhullPoint's printPoint
|
|
|
+*/
|
|
|
+void qh_printpoint(qhT *qh, FILE *fp, const char *string, pointT *point) {
|
|
|
+ int id= qh_pointid(qh, point);
|
|
|
+
|
|
|
+ qh_printpointid(qh, fp, string, qh->hull_dim, point, id);
|
|
|
+} /* printpoint */
|
|
|
+
|
|
|
+void qh_printpointid(qhT *qh, FILE *fp, const char *string, int dim, pointT *point, int id) {
|
|
|
+ int k;
|
|
|
+ realT r; /*bug fix*/
|
|
|
+
|
|
|
+ if (!point)
|
|
|
+ return;
|
|
|
+ if (string) {
|
|
|
+ qh_fprintf(qh, fp, 9211, "%s", string);
|
|
|
+ if (id != qh_IDunknown && id != qh_IDnone)
|
|
|
+ qh_fprintf(qh, fp, 9212, " p%d: ", id);
|
|
|
+ }
|
|
|
+ for (k=dim; k--; ) {
|
|
|
+ r= *point++;
|
|
|
+ if (string)
|
|
|
+ qh_fprintf(qh, fp, 9213, " %8.4g", r);
|
|
|
+ else
|
|
|
+ qh_fprintf(qh, fp, 9214, qh_REAL_1, r);
|
|
|
+ }
|
|
|
+ qh_fprintf(qh, fp, 9215, "\n");
|
|
|
+} /* printpointid */
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printpoint3">-</a>
|
|
|
+
|
|
|
+ qh_printpoint3(qh, fp, point )
|
|
|
+ prints 2-d, 3-d, or 4-d point as Geomview 3-d coordinates
|
|
|
+*/
|
|
|
+void qh_printpoint3(qhT *qh, FILE *fp, pointT *point) {
|
|
|
+ int k;
|
|
|
+ realT p[4];
|
|
|
+
|
|
|
+ qh_projectdim3(qh, point, p);
|
|
|
+ for (k=0; k < 3; k++)
|
|
|
+ qh_fprintf(qh, fp, 9216, "%8.4g ", p[k]);
|
|
|
+ qh_fprintf(qh, fp, 9217, " # p%d\n", qh_pointid(qh, point));
|
|
|
+} /* printpoint3 */
|
|
|
+
|
|
|
+/*----------------------------------------
|
|
|
+-printpoints- print pointids for a set of points starting at index
|
|
|
+ see geom_r.c
|
|
|
+*/
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printpoints_out">-</a>
|
|
|
+
|
|
|
+ qh_printpoints_out(qh, 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(qhT *qh, FILE *fp, facetT *facetlist, setT *facets, boolT printall) {
|
|
|
+ int allpoints= qh->num_points + qh_setsize(qh, 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(qh, allpoints);
|
|
|
+ qh_setzero(qh, points, 0, allpoints);
|
|
|
+ vertices= qh_facetvertices(qh, facetlist, facets, printall);
|
|
|
+ FOREACHvertex_(vertices) {
|
|
|
+ id= qh_pointid(qh, vertex->point);
|
|
|
+ if (id >= 0)
|
|
|
+ SETelem_(points, id)= vertex->point;
|
|
|
+ }
|
|
|
+ if (qh->KEEPinside || qh->KEEPcoplanar || qh->KEEPnearinside) {
|
|
|
+ FORALLfacet_(facetlist) {
|
|
|
+ if (!printall && qh_skipfacet(qh, facet))
|
|
|
+ continue;
|
|
|
+ FOREACHpoint_(facet->coplanarset) {
|
|
|
+ id= qh_pointid(qh, point);
|
|
|
+ if (id >= 0)
|
|
|
+ SETelem_(points, id)= point;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ FOREACHfacet_(facets) {
|
|
|
+ if (!printall && qh_skipfacet(qh, facet))
|
|
|
+ continue;
|
|
|
+ FOREACHpoint_(facet->coplanarset) {
|
|
|
+ id= qh_pointid(qh, point);
|
|
|
+ if (id >= 0)
|
|
|
+ SETelem_(points, id)= point;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ qh_settempfree(qh, &vertices);
|
|
|
+ FOREACHpoint_i_(qh, points) {
|
|
|
+ if (point)
|
|
|
+ numpoints++;
|
|
|
+ }
|
|
|
+ if (qh->CDDoutput)
|
|
|
+ qh_fprintf(qh, fp, 9218, "%s | %s\nbegin\n%d %d real\n", qh->rbox_command,
|
|
|
+ qh->qhull_command, numpoints, qh->hull_dim + 1);
|
|
|
+ else
|
|
|
+ qh_fprintf(qh, fp, 9219, "%d\n%d\n", qh->hull_dim, numpoints);
|
|
|
+ FOREACHpoint_i_(qh, points) {
|
|
|
+ if (point) {
|
|
|
+ if (qh->CDDoutput)
|
|
|
+ qh_fprintf(qh, fp, 9220, "1 ");
|
|
|
+ qh_printpoint(qh, fp, NULL, point);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ if (qh->CDDoutput)
|
|
|
+ qh_fprintf(qh, fp, 9221, "end\n");
|
|
|
+ qh_settempfree(qh, &points);
|
|
|
+} /* printpoints_out */
|
|
|
+
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printpointvect">-</a>
|
|
|
+
|
|
|
+ qh_printpointvect(qh, 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(qhT *qh, 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(qh, diff, qh->hull_dim, True, NULL, NULL);
|
|
|
+ for (k=qh->hull_dim; k--; )
|
|
|
+ pointA[k]= point[k]+diff[k] * radius;
|
|
|
+ qh_printline3geom(qh, fp, point, pointA, color);
|
|
|
+} /* printpointvect */
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printpointvect2">-</a>
|
|
|
+
|
|
|
+ qh_printpointvect2(qh, 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(qhT *qh, FILE *fp, pointT *point, coordT *normal, pointT *center, realT radius) {
|
|
|
+ realT red[3]={1, 0, 0}, yellow[3]={1, 1, 0};
|
|
|
+
|
|
|
+ qh_printpointvect(qh, fp, point, normal, center, radius, red);
|
|
|
+ qh_printpointvect(qh, fp, point, normal, center, -radius, yellow);
|
|
|
+} /* printpointvect2 */
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printridge">-</a>
|
|
|
+
|
|
|
+ qh_printridge(qh, fp, ridge )
|
|
|
+ prints the information in a ridge
|
|
|
+
|
|
|
+ notes:
|
|
|
+ for qh_printfacetridges()
|
|
|
+ same as operator<< [QhullRidge.cpp]
|
|
|
+*/
|
|
|
+void qh_printridge(qhT *qh, FILE *fp, ridgeT *ridge) {
|
|
|
+
|
|
|
+ qh_fprintf(qh, fp, 9222, " - r%d", ridge->id);
|
|
|
+ if (ridge->tested)
|
|
|
+ qh_fprintf(qh, fp, 9223, " tested");
|
|
|
+ if (ridge->nonconvex)
|
|
|
+ qh_fprintf(qh, fp, 9224, " nonconvex");
|
|
|
+ if (ridge->mergevertex)
|
|
|
+ qh_fprintf(qh, fp, 9421, " mergevertex");
|
|
|
+ if (ridge->mergevertex2)
|
|
|
+ qh_fprintf(qh, fp, 9422, " mergevertex2");
|
|
|
+ if (ridge->simplicialtop)
|
|
|
+ qh_fprintf(qh, fp, 9425, " simplicialtop");
|
|
|
+ if (ridge->simplicialbot)
|
|
|
+ qh_fprintf(qh, fp, 9423, " simplicialbot");
|
|
|
+ qh_fprintf(qh, fp, 9225, "\n");
|
|
|
+ qh_printvertices(qh, fp, " vertices:", ridge->vertices);
|
|
|
+ if (ridge->top && ridge->bottom)
|
|
|
+ qh_fprintf(qh, fp, 9226, " between f%d and f%d\n",
|
|
|
+ ridge->top->id, ridge->bottom->id);
|
|
|
+} /* printridge */
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printspheres">-</a>
|
|
|
+
|
|
|
+ qh_printspheres(qh, fp, vertices, radius )
|
|
|
+ prints 3-d vertices as OFF spheres
|
|
|
+
|
|
|
+ notes:
|
|
|
+ inflated octahedron from Stuart Levy earth/mksphere2
|
|
|
+*/
|
|
|
+void qh_printspheres(qhT *qh, FILE *fp, setT *vertices, realT radius) {
|
|
|
+ vertexT *vertex, **vertexp;
|
|
|
+
|
|
|
+ qh->printoutnum++;
|
|
|
+ qh_fprintf(qh, fp, 9227, "{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) {
|
|
|
+ qh_fprintf(qh, fp, 9228, "%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(qh, fp, vertex->point);
|
|
|
+ qh_fprintf(qh, fp, 9229, "1\n");
|
|
|
+ }
|
|
|
+ qh_fprintf(qh, fp, 9230, "}}}\n");
|
|
|
+} /* printspheres */
|
|
|
+
|
|
|
+
|
|
|
+/*----------------------------------------------
|
|
|
+-printsummary-
|
|
|
+ see libqhull_r.c
|
|
|
+*/
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printvdiagram">-</a>
|
|
|
+
|
|
|
+ qh_printvdiagram(qh, 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(qhT *qh, FILE *fp, qh_PRINT 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 {
|
|
|
+ qh_fprintf(qh, qh->ferr, 6219, "qhull internal error (qh_printvdiagram): unknown print format %d.\n", format);
|
|
|
+ qh_errexit(qh, qh_ERRqhull, NULL, NULL);
|
|
|
+ }
|
|
|
+ vertices= qh_markvoronoi(qh, facetlist, facets, printall, &isLower, &numcenters);
|
|
|
+ totcount= qh_printvdiagram2(qh, NULL, NULL, vertices, innerouter, False);
|
|
|
+ qh_fprintf(qh, fp, 9231, "%d\n", totcount);
|
|
|
+ totcount= qh_printvdiagram2(qh, fp, printvridge, vertices, innerouter, True /* inorder*/);
|
|
|
+ qh_settempfree(qh, &vertices);
|
|
|
+#if 0 /* for testing qh_eachvoronoi_all */
|
|
|
+ qh_fprintf(qh, fp, 9232, "\n");
|
|
|
+ totcount= qh_eachvoronoi_all(qh, fp, printvridge, qh->UPPERdelaunay, innerouter, True /* inorder*/);
|
|
|
+ qh_fprintf(qh, fp, 9233, "%d\n", totcount);
|
|
|
+#endif
|
|
|
+} /* printvdiagram */
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printvdiagram2">-</a>
|
|
|
+
|
|
|
+ qh_printvdiagram2(qh, 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(qhT *qh, 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_(qh, vertices) {
|
|
|
+ if (vertex) {
|
|
|
+ if (qh->GOODvertex > 0 && qh_pointid(qh, vertex->point)+1 != qh->GOODvertex)
|
|
|
+ continue;
|
|
|
+ totcount += qh_eachvoronoi(qh, fp, printvridge, vertex, !qh_ALL, innerouter, inorder);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ return totcount;
|
|
|
+} /* printvdiagram2 */
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printvertex">-</a>
|
|
|
+
|
|
|
+ qh_printvertex(qh, fp, vertex )
|
|
|
+ prints the information in a vertex
|
|
|
+ Duplicated as operator<< [QhullVertex.cpp]
|
|
|
+*/
|
|
|
+void qh_printvertex(qhT *qh, FILE *fp, vertexT *vertex) {
|
|
|
+ pointT *point;
|
|
|
+ int k, count= 0;
|
|
|
+ facetT *neighbor, **neighborp;
|
|
|
+ realT r; /*bug fix*/
|
|
|
+
|
|
|
+ if (!vertex) {
|
|
|
+ qh_fprintf(qh, fp, 9234, " NULLvertex\n");
|
|
|
+ return;
|
|
|
+ }
|
|
|
+ qh_fprintf(qh, fp, 9235, "- p%d(v%d):", qh_pointid(qh, vertex->point), vertex->id);
|
|
|
+ point= vertex->point;
|
|
|
+ if (point) {
|
|
|
+ for (k=qh->hull_dim; k--; ) {
|
|
|
+ r= *point++;
|
|
|
+ qh_fprintf(qh, fp, 9236, " %5.2g", r);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ if (vertex->deleted)
|
|
|
+ qh_fprintf(qh, fp, 9237, " deleted");
|
|
|
+ if (vertex->delridge)
|
|
|
+ qh_fprintf(qh, fp, 9238, " delridge");
|
|
|
+ if (vertex->newfacet)
|
|
|
+ qh_fprintf(qh, fp, 9415, " newfacet");
|
|
|
+ if (vertex->seen && qh->IStracing)
|
|
|
+ qh_fprintf(qh, fp, 9416, " seen");
|
|
|
+ if (vertex->seen2 && qh->IStracing)
|
|
|
+ qh_fprintf(qh, fp, 9417, " seen2");
|
|
|
+ qh_fprintf(qh, fp, 9239, "\n");
|
|
|
+ if (vertex->neighbors) {
|
|
|
+ qh_fprintf(qh, fp, 9240, " neighbors:");
|
|
|
+ FOREACHneighbor_(vertex) {
|
|
|
+ if (++count % 100 == 0)
|
|
|
+ qh_fprintf(qh, fp, 9241, "\n ");
|
|
|
+ qh_fprintf(qh, fp, 9242, " f%d", neighbor->id);
|
|
|
+ }
|
|
|
+ qh_fprintf(qh, fp, 9243, "\n");
|
|
|
+ }
|
|
|
+} /* printvertex */
|
|
|
+
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printvertexlist">-</a>
|
|
|
+
|
|
|
+ qh_printvertexlist(qh, fp, string, facetlist, facets, printall )
|
|
|
+ prints vertices used by a facetlist or facet set
|
|
|
+ tests qh_skipfacet() if !printall
|
|
|
+*/
|
|
|
+void qh_printvertexlist(qhT *qh, FILE *fp, const char* string, facetT *facetlist,
|
|
|
+ setT *facets, boolT printall) {
|
|
|
+ vertexT *vertex, **vertexp;
|
|
|
+ setT *vertices;
|
|
|
+
|
|
|
+ vertices= qh_facetvertices(qh, facetlist, facets, printall);
|
|
|
+ qh_fprintf(qh, fp, 9244, "%s", string);
|
|
|
+ FOREACHvertex_(vertices)
|
|
|
+ qh_printvertex(qh, fp, vertex);
|
|
|
+ qh_settempfree(qh, &vertices);
|
|
|
+} /* printvertexlist */
|
|
|
+
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printvertices">-</a>
|
|
|
+
|
|
|
+ qh_printvertices(qh, fp, string, vertices )
|
|
|
+ prints vertices in a set
|
|
|
+ duplicated as printVertexSet [QhullVertex.cpp]
|
|
|
+*/
|
|
|
+void qh_printvertices(qhT *qh, FILE *fp, const char* string, setT *vertices) {
|
|
|
+ vertexT *vertex, **vertexp;
|
|
|
+
|
|
|
+ qh_fprintf(qh, fp, 9245, "%s", string);
|
|
|
+ FOREACHvertex_(vertices)
|
|
|
+ qh_fprintf(qh, fp, 9246, " p%d(v%d)", qh_pointid(qh, vertex->point), vertex->id);
|
|
|
+ qh_fprintf(qh, fp, 9247, "\n");
|
|
|
+} /* printvertices */
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printvneighbors">-</a>
|
|
|
+
|
|
|
+ qh_printvneighbors(qh, 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(qhT *qh, 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, qh->other_points);
|
|
|
+ vertexT *vertex, **vertexp;
|
|
|
+ int vertex_i, vertex_n;
|
|
|
+ facetT *facet, **facetp, *neighbor, **neighborp;
|
|
|
+ pointT *point, **pointp;
|
|
|
+
|
|
|
+ qh_countfacets(qh, facetlist, facets, printall, &numfacets, &numsimplicial,
|
|
|
+ &totneighbors, &numridges, &numcoplanars, &numtricoplanars); /* sets facet->visitid */
|
|
|
+ qh_fprintf(qh, fp, 9248, "%d\n", numpoints);
|
|
|
+ qh_vertexneighbors(qh);
|
|
|
+ vertices= qh_facetvertices(qh, facetlist, facets, printall);
|
|
|
+ vertex_points= qh_settemp(qh, numpoints);
|
|
|
+ coplanar_points= qh_settemp(qh, numpoints);
|
|
|
+ qh_setzero(qh, vertex_points, 0, numpoints);
|
|
|
+ qh_setzero(qh, coplanar_points, 0, numpoints);
|
|
|
+ FOREACHvertex_(vertices)
|
|
|
+ qh_point_add(qh, vertex_points, vertex->point, vertex);
|
|
|
+ FORALLfacet_(facetlist) {
|
|
|
+ FOREACHpoint_(facet->coplanarset)
|
|
|
+ qh_point_add(qh, coplanar_points, point, facet);
|
|
|
+ }
|
|
|
+ FOREACHfacet_(facets) {
|
|
|
+ FOREACHpoint_(facet->coplanarset)
|
|
|
+ qh_point_add(qh, coplanar_points, point, facet);
|
|
|
+ }
|
|
|
+ FOREACHvertex_i_(qh, vertex_points) {
|
|
|
+ if (vertex) {
|
|
|
+ numneighbors= qh_setsize(qh, vertex->neighbors);
|
|
|
+ qh_fprintf(qh, fp, 9249, "%d", numneighbors);
|
|
|
+ qh_order_vertexneighbors(qh, vertex);
|
|
|
+ FOREACHneighbor_(vertex)
|
|
|
+ qh_fprintf(qh, fp, 9250, " %d",
|
|
|
+ neighbor->visitid ? neighbor->visitid - 1 : 0 - neighbor->id);
|
|
|
+ qh_fprintf(qh, fp, 9251, "\n");
|
|
|
+ }else if ((facet= SETelemt_(coplanar_points, vertex_i, facetT)))
|
|
|
+ qh_fprintf(qh, fp, 9252, "1 %d\n",
|
|
|
+ facet->visitid ? facet->visitid - 1 : 0 - facet->id);
|
|
|
+ else
|
|
|
+ qh_fprintf(qh, fp, 9253, "0\n");
|
|
|
+ }
|
|
|
+ qh_settempfree(qh, &coplanar_points);
|
|
|
+ qh_settempfree(qh, &vertex_points);
|
|
|
+ qh_settempfree(qh, &vertices);
|
|
|
+} /* printvneighbors */
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printvoronoi">-</a>
|
|
|
+
|
|
|
+ qh_printvoronoi(qh, 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(qhT *qh, FILE *fp, qh_PRINT 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(qh, facetlist, facets, printall, &isLower, &numcenters);
|
|
|
+ FOREACHvertex_i_(qh, 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)
|
|
|
+ qh_fprintf(qh, fp, 9254, "{appearance {+edge -face} OFF %d %d 1 # Voronoi centers and cells\n",
|
|
|
+ numcenters, numvertices);
|
|
|
+ else
|
|
|
+ qh_fprintf(qh, fp, 9255, "%d\n%d %d 1\n", qh->hull_dim-1, numcenters, qh_setsize(qh, vertices));
|
|
|
+ if (format == qh_PRINTgeom) {
|
|
|
+ for (k=qh->hull_dim-1; k--; )
|
|
|
+ qh_fprintf(qh, fp, 9256, qh_REAL_1, 0.0);
|
|
|
+ qh_fprintf(qh, fp, 9257, " 0 # infinity not used\n");
|
|
|
+ }else {
|
|
|
+ for (k=qh->hull_dim-1; k--; )
|
|
|
+ qh_fprintf(qh, fp, 9258, qh_REAL_1, qh_INFINITE);
|
|
|
+ qh_fprintf(qh, fp, 9259, "\n");
|
|
|
+ }
|
|
|
+ FORALLfacet_(facetlist) {
|
|
|
+ if (facet->visitid && facet->visitid < numfacets) {
|
|
|
+ if (format == qh_PRINTgeom)
|
|
|
+ qh_fprintf(qh, fp, 9260, "# %d f%d\n", vid++, facet->id);
|
|
|
+ qh_printcenter(qh, fp, format, NULL, facet);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ FOREACHfacet_(facets) {
|
|
|
+ if (facet->visitid && facet->visitid < numfacets) {
|
|
|
+ if (format == qh_PRINTgeom)
|
|
|
+ qh_fprintf(qh, fp, 9261, "# %d f%d\n", vid++, facet->id);
|
|
|
+ qh_printcenter(qh, fp, format, NULL, facet);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ FOREACHvertex_i_(qh, vertices) {
|
|
|
+ numneighbors= 0;
|
|
|
+ numinf=0;
|
|
|
+ if (vertex) {
|
|
|
+ qh_order_vertexneighbors(qh, vertex);
|
|
|
+ FOREACHneighbor_(vertex) {
|
|
|
+ if (neighbor->visitid == 0)
|
|
|
+ numinf= 1;
|
|
|
+ else if (neighbor->visitid < numfacets)
|
|
|
+ numneighbors++;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ if (format == qh_PRINTgeom) {
|
|
|
+ if (vertex) {
|
|
|
+ qh_fprintf(qh, fp, 9262, "%d", numneighbors);
|
|
|
+ FOREACHneighbor_(vertex) {
|
|
|
+ if (neighbor->visitid && neighbor->visitid < numfacets)
|
|
|
+ qh_fprintf(qh, fp, 9263, " %d", neighbor->visitid);
|
|
|
+ }
|
|
|
+ qh_fprintf(qh, fp, 9264, " # p%d(v%d)\n", vertex_i, vertex->id);
|
|
|
+ }else
|
|
|
+ qh_fprintf(qh, fp, 9265, " # p%d is coplanar or isolated\n", vertex_i);
|
|
|
+ }else {
|
|
|
+ if (numinf)
|
|
|
+ numneighbors++;
|
|
|
+ qh_fprintf(qh, fp, 9266, "%d", numneighbors);
|
|
|
+ if (vertex) {
|
|
|
+ FOREACHneighbor_(vertex) {
|
|
|
+ if (neighbor->visitid == 0) {
|
|
|
+ if (numinf) {
|
|
|
+ numinf= 0;
|
|
|
+ qh_fprintf(qh, fp, 9267, " %d", neighbor->visitid);
|
|
|
+ }
|
|
|
+ }else if (neighbor->visitid < numfacets)
|
|
|
+ qh_fprintf(qh, fp, 9268, " %d", neighbor->visitid);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ qh_fprintf(qh, fp, 9269, "\n");
|
|
|
+ }
|
|
|
+ }
|
|
|
+ if (format == qh_PRINTgeom)
|
|
|
+ qh_fprintf(qh, fp, 9270, "}\n");
|
|
|
+ qh_settempfree(qh, &vertices);
|
|
|
+} /* printvoronoi */
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printvnorm">-</a>
|
|
|
+
|
|
|
+ qh_printvnorm(qh, 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
|
|
|
+
|
|
|
+ note:
|
|
|
+ parameter unbounded is UNUSED by this callback
|
|
|
+
|
|
|
+ see:
|
|
|
+ qh_printvdiagram()
|
|
|
+ qh_eachvoronoi()
|
|
|
+*/
|
|
|
+void qh_printvnorm(qhT *qh, FILE *fp, vertexT *vertex, vertexT *vertexA, setT *centers, boolT unbounded) {
|
|
|
+ pointT *normal;
|
|
|
+ realT offset;
|
|
|
+ int k;
|
|
|
+ QHULL_UNUSED(unbounded);
|
|
|
+
|
|
|
+ normal= qh_detvnorm(qh, vertex, vertexA, centers, &offset);
|
|
|
+ qh_fprintf(qh, fp, 9271, "%d %d %d ",
|
|
|
+ 2+qh->hull_dim, qh_pointid(qh, vertex->point), qh_pointid(qh, vertexA->point));
|
|
|
+ for (k=0; k< qh->hull_dim-1; k++)
|
|
|
+ qh_fprintf(qh, fp, 9272, qh_REAL_1, normal[k]);
|
|
|
+ qh_fprintf(qh, fp, 9273, qh_REAL_1, offset);
|
|
|
+ qh_fprintf(qh, fp, 9274, "\n");
|
|
|
+} /* printvnorm */
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="printvridge">-</a>
|
|
|
+
|
|
|
+ qh_printvridge(qh, 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
|
|
|
+ parameter unbounded is UNUSED
|
|
|
+*/
|
|
|
+void qh_printvridge(qhT *qh, FILE *fp, vertexT *vertex, vertexT *vertexA, setT *centers, boolT unbounded) {
|
|
|
+ facetT *facet, **facetp;
|
|
|
+ QHULL_UNUSED(unbounded);
|
|
|
+
|
|
|
+ qh_fprintf(qh, fp, 9275, "%d %d %d", qh_setsize(qh, centers)+2,
|
|
|
+ qh_pointid(qh, vertex->point), qh_pointid(qh, vertexA->point));
|
|
|
+ FOREACHfacet_(centers)
|
|
|
+ qh_fprintf(qh, fp, 9276, " %d", facet->visitid);
|
|
|
+ qh_fprintf(qh, fp, 9277, "\n");
|
|
|
+} /* printvridge */
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="projectdim3">-</a>
|
|
|
+
|
|
|
+ qh_projectdim3(qh, 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(qhT *qh, 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_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="readfeasible">-</a>
|
|
|
+
|
|
|
+ qh_readfeasible(qh, dim, curline )
|
|
|
+ read feasible point from current line and qh.fin
|
|
|
+
|
|
|
+ returns:
|
|
|
+ number of lines read from qh.fin
|
|
|
+ sets qh.feasible_point with malloc'd coordinates
|
|
|
+
|
|
|
+ notes:
|
|
|
+ checks for qh.HALFspace
|
|
|
+ assumes dim > 1
|
|
|
+
|
|
|
+ see:
|
|
|
+ qh_setfeasible
|
|
|
+*/
|
|
|
+int qh_readfeasible(qhT *qh, int dim, const char *curline) {
|
|
|
+ boolT isfirst= True;
|
|
|
+ int linecount= 0, tokcount= 0;
|
|
|
+ const char *s;
|
|
|
+ char *t, firstline[qh_MAXfirst+1];
|
|
|
+ coordT *coords, value;
|
|
|
+
|
|
|
+ if (!qh->HALFspace) {
|
|
|
+ qh_fprintf(qh, qh->ferr, 6070, "qhull input error: feasible point(dim 1 coords) is only valid for halfspace intersection\n");
|
|
|
+ qh_errexit(qh, qh_ERRinput, NULL, NULL);
|
|
|
+ }
|
|
|
+ if (qh->feasible_string)
|
|
|
+ qh_fprintf(qh, qh->ferr, 7057, "qhull input warning: feasible point(dim 1 coords) overrides 'Hn,n,n' feasible point for halfspace intersection\n");
|
|
|
+ if (!(qh->feasible_point= (coordT *)qh_malloc((size_t)dim * sizeof(coordT)))) {
|
|
|
+ qh_fprintf(qh, qh->ferr, 6071, "qhull error: insufficient memory for feasible point\n");
|
|
|
+ qh_errexit(qh, qh_ERRmem, NULL, NULL);
|
|
|
+ }
|
|
|
+ coords= qh->feasible_point;
|
|
|
+ while ((s= (isfirst ? curline : 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) {
|
|
|
+ qh_fprintf(qh, qh->ferr, 6072, "qhull input error: coordinates for feasible point do not finish out the line: %s\n",
|
|
|
+ s);
|
|
|
+ qh_errexit(qh, qh_ERRinput, NULL, NULL);
|
|
|
+ }
|
|
|
+ return linecount;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ qh_fprintf(qh, qh->ferr, 6073, "qhull input error: only %d coordinates. Could not read %d-d feasible point.\n",
|
|
|
+ tokcount, dim);
|
|
|
+ qh_errexit(qh, qh_ERRinput, NULL, NULL);
|
|
|
+ return 0;
|
|
|
+} /* readfeasible */
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="readpoints">-</a>
|
|
|
+
|
|
|
+ qh_readpoints(qh, 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
|
|
|
+ QH11012 FIX: qh_readpoints needs rewriting, too long
|
|
|
+*/
|
|
|
+coordT *qh_readpoints(qhT *qh, 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= 0, *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(qh, dimfeasible, s);
|
|
|
+ else
|
|
|
+ dimfeasible= 0;
|
|
|
+ }else if (!memcmp(firstline, "begin", (size_t)5) || !memcmp(firstline, "BEGIN", (size_t)5))
|
|
|
+ break;
|
|
|
+ else if (!*qh->rbox_command)
|
|
|
+ strncat(qh->rbox_command, s, sizeof(qh->rbox_command)-1);
|
|
|
+ }
|
|
|
+ if (!s) {
|
|
|
+ qh_fprintf(qh, qh->ferr, 6074, "qhull input error: missing \"begin\" for cdd-formated input\n");
|
|
|
+ qh_errexit(qh, qh_ERRinput, NULL, NULL);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ while (!numinput && (s= fgets(firstline, qh_MAXfirst, qh->fin))) {
|
|
|
+ linecount++;
|
|
|
+ if (!memcmp(s, "begin", (size_t)5) || !memcmp(s, "BEGIN", (size_t)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(qh, diminput, s); /* checks if ok */
|
|
|
+ dimfeasible= diminput;
|
|
|
+ diminput= numinput= 0;
|
|
|
+ }else
|
|
|
+ break;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ if (!s) {
|
|
|
+ qh_fprintf(qh, qh->ferr, 6075, "qhull input error: short input file. Did not find dimension and number of points\n");
|
|
|
+ qh_errexit(qh, 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) {
|
|
|
+ qh_fprintf(qh, qh->ferr, 6220, "qhull input error: dimension %d (first or smaller number) should be at least 2\n",
|
|
|
+ diminput);
|
|
|
+ qh_errexit(qh, qh_ERRinput, NULL, NULL);
|
|
|
+ }
|
|
|
+ if (numinput < 1 || numinput > qh_POINTSmax) {
|
|
|
+ qh_fprintf(qh, qh->ferr, 6411, "qhull input error: expecting between 1 and %d points. Got %d %d-d points\n",
|
|
|
+ qh_POINTSmax, numinput, diminput);
|
|
|
+ qh_errexit(qh, qh_ERRinput, NULL, NULL);
|
|
|
+ /* same error message in qh_initqhull_globals */
|
|
|
+ }
|
|
|
+
|
|
|
+ if (isdelaunay && qh->HALFspace) {
|
|
|
+ qh_fprintf(qh, qh->ferr, 6037, "qhull option error (qh_readpoints): can not use Delaunay('d') or Voronoi('v') with halfspace intersection('H')\n");
|
|
|
+ qh_errexit(qh, qh_ERRinput, NULL, NULL);
|
|
|
+ /* otherwise corrupted memory allocations, same error message as in qh_initqhull_globals */
|
|
|
+ }else 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) {
|
|
|
+ qh_fprintf(qh, qh->ferr, 6221, "qhull input error: dimension %d (first number, includes offset) should be at least 3 for halfspaces\n",
|
|
|
+ diminput);
|
|
|
+ qh_errexit(qh, qh_ERRinput, NULL, NULL);
|
|
|
+ }
|
|
|
+ if (dimfeasible) {
|
|
|
+ if (dimfeasible != *dimension) {
|
|
|
+ qh_fprintf(qh, qh->ferr, 6222, "qhull input error: dimension %d of feasible point is not one less than dimension %d for halfspaces\n",
|
|
|
+ dimfeasible, diminput);
|
|
|
+ qh_errexit(qh, qh_ERRinput, NULL, NULL);
|
|
|
+ }
|
|
|
+ }else
|
|
|
+ qh_setfeasible(qh, *dimension);
|
|
|
+ }else {
|
|
|
+ if (qh->CDDinput)
|
|
|
+ *dimension= diminput-1;
|
|
|
+ else
|
|
|
+ *dimension= diminput;
|
|
|
+ *numpoints= numinput;
|
|
|
+ }
|
|
|
+ qh->normal_size= *dimension * (int)sizeof(coordT); /* for tracing with qh_printpoint */
|
|
|
+ if (qh->HALFspace) {
|
|
|
+ qh->half_space= coordp= (coordT *)qh_malloc((size_t)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 *)qh_malloc((size_t)(qh->maxline+1) * sizeof(char));
|
|
|
+ *ismalloc= True; /* use malloc since memory not setup */
|
|
|
+ coords= points= qh->temp_malloc= /* numinput and diminput >=2 by QH6220 */
|
|
|
+ (coordT *)qh_malloc((size_t)((*numpoints)*(*dimension))*sizeof(coordT));
|
|
|
+ if (!coords || !qh->line || (qh->HALFspace && !qh->half_space)) {
|
|
|
+ qh_fprintf(qh, qh->ferr, 6076, "qhull error: insufficient memory to read %d points\n",
|
|
|
+ numinput);
|
|
|
+ qh_errexit(qh, 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", (size_t)3) || !memcmp(s, "END", (size_t)3)) {
|
|
|
+ if (qh->CDDinput )
|
|
|
+ break;
|
|
|
+ else if (wasbegin)
|
|
|
+ qh_fprintf(qh, qh->ferr, 7058, "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) {
|
|
|
+ qh_fprintf(qh, qh->ferr, 6077, "qhull input error: for cdd format, point at line %d does not start with '1'\n",
|
|
|
+ linecount);
|
|
|
+ qh_errexit(qh, 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(qh, *dimension, coords, &coords, normalp, offsetp, qh->feasible_point)) {
|
|
|
+ qh_fprintf(qh, qh->ferr, 8048, "The halfspace was on line %d\n", linecount);
|
|
|
+ if (wasbegin)
|
|
|
+ qh_fprintf(qh, qh->ferr, 8049, "The input appears to be in cdd format. If so, you should use option 'Fd'\n");
|
|
|
+ qh_errexit(qh, 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) {
|
|
|
+ qh_fprintf(qh, qh->ferr, 6078, "qhull input error: line %d contained more than %d characters\n",
|
|
|
+ linecount, (int) (s - qh->line)); /* WARN64 */
|
|
|
+ qh_errexit(qh, qh_ERRinput, NULL, NULL);
|
|
|
+ }
|
|
|
+ isfirst= False;
|
|
|
+ }
|
|
|
+ if (qh->rbox_command[0])
|
|
|
+ qh->rbox_command[strlen(qh->rbox_command)-1]= '\0'; /* remove \n, previous qh_errexit's display command as two lines */
|
|
|
+ if (tokcount != maxcount) {
|
|
|
+ newnum= fmin_(numinput, tokcount/diminput);
|
|
|
+ if (qh->ALLOWshort)
|
|
|
+ qh_fprintf(qh, qh->ferr, 7073, "qhull warning: instead of %d points in %d-d, input contains %d points and %d extra coordinates.\n",
|
|
|
+ numinput, diminput, tokcount/diminput, tokcount % diminput);
|
|
|
+ else
|
|
|
+ qh_fprintf(qh, qh->ferr, 6410, "qhull error: instead of %d points in %d-d, input contains %d points and %d extra coordinates.\n",
|
|
|
+ numinput, diminput, tokcount/diminput, tokcount % diminput);
|
|
|
+ if (firsttext)
|
|
|
+ qh_fprintf(qh, qh->ferr, 8051, " Line %d is the first comment.\n", firsttext);
|
|
|
+ qh_fprintf(qh, qh->ferr, 8033, " Line %d is the first point.\n", firstpoint);
|
|
|
+ if (firstshort)
|
|
|
+ qh_fprintf(qh, qh->ferr, 8052, " Line %d is the first short line.\n", firstshort);
|
|
|
+ if (firstlong)
|
|
|
+ qh_fprintf(qh, qh->ferr, 8053, " Line %d is the first long line.\n", firstlong);
|
|
|
+ if (qh->ALLOWshort)
|
|
|
+ qh_fprintf(qh, qh->ferr, 8054, " Continuing with %d points.\n", newnum);
|
|
|
+ else {
|
|
|
+ qh_fprintf(qh, qh->ferr, 8077, " Override with option 'Qa' (allow-short)\n");
|
|
|
+ qh_errexit(qh, qh_ERRinput, NULL, NULL);
|
|
|
+ }
|
|
|
+ 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 (!strcmp(qh->rbox_command, "./rbox D4"))
|
|
|
+ qh_fprintf(qh, qh->ferr, 8055, "\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_r.h. If you discover the source of the problem,\n\
|
|
|
+please send mail to qhull_bug@qhull.org.\n\
|
|
|
+\n\
|
|
|
+Type 'qhull' for a short list of options.\n");
|
|
|
+ qh_free(qh->line);
|
|
|
+ qh->line= NULL;
|
|
|
+ if (qh->half_space) {
|
|
|
+ qh_free(qh->half_space);
|
|
|
+ qh->half_space= NULL;
|
|
|
+ }
|
|
|
+ qh->temp_malloc= NULL;
|
|
|
+ trace1((qh, qh->ferr, 1008,"qh_readpoints: read in %d %d-dimensional points\n",
|
|
|
+ numinput, diminput));
|
|
|
+ return(points);
|
|
|
+} /* readpoints */
|
|
|
+
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="setfeasible">-</a>
|
|
|
+
|
|
|
+ qh_setfeasible(qh, dim )
|
|
|
+ set qh.feasible_point 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()
|
|
|
+ called only once from qh_new_qhull, otherwise leaks memory
|
|
|
+*/
|
|
|
+void qh_setfeasible(qhT *qh, int dim) {
|
|
|
+ int tokcount= 0;
|
|
|
+ char *s;
|
|
|
+ coordT *coords, value;
|
|
|
+
|
|
|
+ if (!(s= qh->feasible_string)) {
|
|
|
+ qh_fprintf(qh, qh->ferr, 6223, "qhull input error: halfspace intersection needs a feasible point. Either prepend the input with 1 point or use 'Hn,n,n'. See manual.\n");
|
|
|
+ qh_errexit(qh, qh_ERRinput, NULL, NULL);
|
|
|
+ }
|
|
|
+ if (!(qh->feasible_point= (pointT *)qh_malloc((size_t)dim * sizeof(coordT)))) {
|
|
|
+ qh_fprintf(qh, qh->ferr, 6079, "qhull error: insufficient memory for 'Hn,n,n'\n");
|
|
|
+ qh_errexit(qh, qh_ERRmem, NULL, NULL);
|
|
|
+ }
|
|
|
+ coords= qh->feasible_point;
|
|
|
+ while (*s) {
|
|
|
+ value= qh_strtod(s, &s);
|
|
|
+ if (++tokcount > dim) {
|
|
|
+ qh_fprintf(qh, qh->ferr, 7059, "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_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="skipfacet">-</a>
|
|
|
+
|
|
|
+ qh_skipfacet(qh, 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(qhT *qh, 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(qh, facet->normal, NULL));
|
|
|
+} /* skipfacet */
|
|
|
+
|
|
|
+/*-<a href="qh-io_r.htm#TOC"
|
|
|
+ >-------------------------------</a><a name="skipfilename">-</a>
|
|
|
+
|
|
|
+ qh_skipfilename(qh, string )
|
|
|
+ returns pointer to character after filename
|
|
|
+
|
|
|
+ notes:
|
|
|
+ skips leading spaces
|
|
|
+ ends with spacing or eol
|
|
|
+ if starts with ' or " ends with the same, skipping \' or \"
|
|
|
+ For qhull, qh_argv_to_command() only uses double quotes
|
|
|
+*/
|
|
|
+char *qh_skipfilename(qhT *qh, char *filename) {
|
|
|
+ char *s= filename; /* non-const due to return */
|
|
|
+ char c;
|
|
|
+
|
|
|
+ while (*s && isspace(*s))
|
|
|
+ s++;
|
|
|
+ c= *s++;
|
|
|
+ if (c == '\0') {
|
|
|
+ qh_fprintf(qh, qh->ferr, 6204, "qhull input error: filename expected, none found.\n");
|
|
|
+ qh_errexit(qh, qh_ERRinput, NULL, NULL);
|
|
|
+ }
|
|
|
+ if (c == '\'' || c == '"') {
|
|
|
+ while (*s !=c || s[-1] == '\\') {
|
|
|
+ if (!*s) {
|
|
|
+ qh_fprintf(qh, qh->ferr, 6203, "qhull input error: missing quote after filename -- %s\n", filename);
|
|
|
+ qh_errexit(qh, qh_ERRinput, NULL, NULL);
|
|
|
+ }
|
|
|
+ s++;
|
|
|
+ }
|
|
|
+ s++;
|
|
|
+ }
|
|
|
+ else while (*s && !isspace(*s))
|
|
|
+ s++;
|
|
|
+ return s;
|
|
|
+} /* skipfilename */
|
|
|
+
|
yuryalekseev