📄 io.c
字号:
/*<html><pre> -<a href="qh-c.htm#io"
>-------------------------------</a><a name="TOP">-</a>
io.c
Input/Output routines of qhull application
see qh-c.htm and io.h
see user.c for qh_errprint and qh_printfacetlist
unix.c calls qh_readpoints and qh_produce_output
unix.c and user.c are the only callers of io.c functions
This allows the user to avoid loading io.o from qhull.a
copyright (c) 1993-1999 The Geometry Center
*/
#include "qhull_a.h"
/*========= -prototypes for internal functions ========= */
static int qh_compare_facetarea(const void *p1, const void *p2);
static int qh_compare_facetmerge(const void *p1, const void *p2);
static int qh_compare_facetvisit(const void *p1, const void *p2);
int qh_compare_vertexpoint(const void *p1, const void *p2); /* not used */
/*========= -functions in alphabetical order after qh_produce_output() =====*/
/*-<a href="qh-c.htm#io"
>-------------------------------</a><a name="produce_output">-</a>
qh_produce_output()
prints out the result of qhull in desired format
if qh.GETarea
computes and prints area and volume
qh.PRINTout[] is an array of output formats
notes:
prints output in qh.PRINTout order
*/
void qh_produce_output(void) {
int i, tempsize= qh_setsize ((setT*)qhmem.tempstack), d_1;
if (qh VORONOI) {
qh_clearcenters (qh_ASvoronoi);
qh_vertexneighbors();
}
if (qh GETarea)
qh_getarea(qh facet_list);
qh_findgood_all (qh facet_list);
if (qh KEEParea || qh KEEPmerge || qh KEEPminArea < REALmax/2)
qh_markkeep (qh facet_list);
if (qh PRINTsummary)
qh_printsummary(qh ferr);
else if (qh PRINTout[0] == qh_PRINTnone)
qh_printsummary(qh fout);
for (i= 0; i < qh_PRINTEND; i++)
qh_printfacets (qh fout, qh PRINTout[i], qh facet_list, NULL, !qh_ALL);
qh_allstatistics();
if (qh PRINTprecision && !qh MERGING && (qh JOGGLEmax > REALmax/2 || qh RERUN))
qh_printstats (qh ferr, qhstat precision, NULL);
if (qh VERIFYoutput && (zzval_(Zridge) > 0 || zzval_(Zridgemid) > 0))
qh_printstats (qh ferr, qhstat vridges, NULL);
if (qh PRINTstatistics) {
qh_collectstatistics();
qh_printstatistics(qh ferr, "");
qh_memstatistics (qh ferr);
d_1= sizeof(setT) + (qh hull_dim - 1) * SETelemsize;
fprintf(qh ferr, "\
size in bytes: hashentry %d merge %d ridge %d vertex %d facet %d\n\
normal %d ridge vertices %d facet vertices or neighbors %d\n",
sizeof(hashentryT), sizeof(mergeT), sizeof(ridgeT),
sizeof(vertexT), sizeof(facetT),
qh normal_size, d_1, d_1 + SETelemsize);
}
if (qh_setsize ((setT*)qhmem.tempstack) != tempsize) {
fprintf (qh ferr, "qhull internal error (qh_produce_output): temporary sets not empty (%d)\n",
qh_setsize ((setT*)qhmem.tempstack));
qh_errexit (qh_ERRqhull, NULL, NULL);
}
} /* produce_output */
/*-<a href="qh-c.htm#io"
>-------------------------------</a><a name="dfacet">-</a>
dfacet( id )
print facet by id, for debugging
*/
void dfacet (unsigned id) {
facetT *facet;
FORALLfacets {
if (facet->id == id) {
qh_printfacet (qh fout, facet);
break;
}
}
} /* dfacet */
/*-<a href="qh-c.htm#io"
>-------------------------------</a><a name="dvertex">-</a>
dvertex( id )
print vertex by id, for debugging
*/
void dvertex (unsigned id) {
vertexT *vertex;
FORALLvertices {
if (vertex->id == id) {
qh_printvertex (qh fout, vertex);
break;
}
}
} /* dvertex */
/*-<a href="qh-c.htm#io"
>-------------------------------</a><a name="compare_vertexpoint">-</a>
qh_compare_vertexpoint( p1, p2 )
used by qsort() to order vertices by point id
*/
int qh_compare_vertexpoint(const void *p1, const void *p2) {
vertexT *a= *((vertexT **)p1), *b= *((vertexT **)p2);
return ((qh_pointid(a->point) > qh_pointid(b->point)?1:-1));
} /* compare_vertexpoint */
/*-<a href="qh-c.htm#io"
>-------------------------------</a><a name="compare_facetarea">-</a>
qh_compare_facetarea( p1, p2 )
used by qsort() to order facets by area
*/
static int qh_compare_facetarea(const void *p1, const void *p2) {
facetT *a= *((facetT **)p1), *b= *((facetT **)p2);
if (!a->isarea)
return -1;
if (!b->isarea)
return 1;
if (a->f.area > b->f.area)
return 1;
else if (a->f.area == b->f.area)
return 0;
return -1;
} /* compare_facetarea */
/*-<a href="qh-c.htm#io"
>-------------------------------</a><a name="compare_facetmerge">-</a>
qh_compare_facetmerge( p1, p2 )
used by qsort() to order facets by number of merges
*/
static int qh_compare_facetmerge(const void *p1, const void *p2) {
facetT *a= *((facetT **)p1), *b= *((facetT **)p2);
return (a->nummerge - b->nummerge);
} /* compare_facetvisit */
/*-<a href="qh-c.htm#io"
>-------------------------------</a><a name="compare_facetvisit">-</a>
qh_compare_facetvisit( p1, p2 )
used by qsort() to order facets by visit id or id
*/
static int qh_compare_facetvisit(const void *p1, const void *p2) {
facetT *a= *((facetT **)p1), *b= *((facetT **)p2);
int i,j;
if (!(i= a->visitid))
i= - a->id; /* do not convert to int */
if (!(j= b->visitid))
j= - b->id;
return (i - j);
} /* compare_facetvisit */
/*-<a href="qh-c.htm#io"
>-------------------------------</a><a name="countfacets">-</a>
qh_countfacets( facetlist, facets, printall,
numfacets, numsimplicial, totneighbors, numridges, numcoplanar )
count good facets for printing and set visitid
if allfacets, ignores qh_skipfacet()
returns:
numfacets, numsimplicial, total neighbors, numridges, coplanars
each facet with ->visitid indicating 1-relative position
->visitid==0 indicates not good
notes
if qh.NEWfacets,
does not count visible facets (matches qh_printafacet)
design:
for all facets on facetlist and in facets set
unless facet is skipped or visible (i.e., will be deleted)
mark facet->visitid
update counts
*/
void qh_countfacets (facetT *facetlist, setT *facets, boolT printall,
int *numfacetsp, int *numsimplicialp, int *totneighborsp, int *numridgesp, int *numcoplanarsp) {
facetT *facet, **facetp;
int numfacets= 0, numsimplicial= 0, numridges= 0, totneighbors= 0, numcoplanars= 0;
FORALLfacet_(facetlist) {
if ((facet->visible && qh NEWfacets)
|| (!printall && qh_skipfacet(facet)))
facet->visitid= 0;
else {
facet->visitid= ++numfacets;
totneighbors += qh_setsize (facet->neighbors);
if (facet->simplicial)
numsimplicial++;
else
numridges += qh_setsize (facet->ridges);
if (facet->coplanarset)
numcoplanars += qh_setsize (facet->coplanarset);
}
}
FOREACHfacet_(facets) {
if ((facet->visible && qh NEWfacets)
|| (!printall && qh_skipfacet(facet)))
facet->visitid= 0;
else {
facet->visitid= ++numfacets;
totneighbors += qh_setsize (facet->neighbors);
if (facet->simplicial)
numsimplicial++;
else
numridges += qh_setsize (facet->ridges);
if (facet->coplanarset)
numcoplanars += qh_setsize (facet->coplanarset);
}
}
qh visit_id += numfacets+1;
*numfacetsp= numfacets;
*numsimplicialp= numsimplicial;
*totneighborsp= totneighbors;
*numridgesp= numridges;
*numcoplanarsp= numcoplanars;
} /* countfacets */
/*-<a href="qh-c.htm#io"
>-------------------------------</a><a name="detvnorm">-</a>
qh_detvnorm( vertex, vertexA, centers, offset )
compute separating plane of the Voronoi diagram for a pair of input sites
centers= set of facets (i.e., Voronoi vertices)
facet->visitid= 0 iff vertex-at-infinity (i.e., unbounded)
assumes:
qh_ASvoronoi and qh_vertexneighbors() already set
returns:
norm
a pointer into qh.gm_matrix to qh.hull_dim-1 reals
copy the data before reusing qh.gm_matrix
offset
if 'QVn'
sign adjusted so that qh.GOODvertexp is inside
else
sign adjusted so that vertex is inside
qh.gm_matrix= simplex of points from centers relative to first center
notes:
in io.c so that code for 'v Tv' can be removed by removing io.c
returns pointer into qh.gm_matrix to avoid tracking of temporary memory
design:
determine midpoint of input sites
build points as the set of Voronoi vertices
select a simplex from points (if necessary)
include midpoint if the Voronoi region is unbounded
relocate the first vertex of the simplex to the origin
compute the normalized hyperplane through the simplex
orient the hyperplane toward 'QVn' or 'vertex'
if 'Tv' or 'Ts'
if bounded
test that hyperplane is the perpendicular bisector of the input sites
test that Voronoi vertices not in the simplex are still on the hyperplane
free up temporary memory
*/
pointT *qh_detvnorm (vertexT *vertex, vertexT *vertexA, setT *centers, realT *offsetp) {
facetT *facet, **facetp;
int i, k, pointid, pointidA, point_i, point_n;
setT *simplex= NULL;
pointT *point, **pointp, *point0, *midpoint, *normal, *inpoint;
coordT *coord, *gmcoord, *normalp;
setT *points= qh_settemp (qh TEMPsize);
boolT nearzero= False;
boolT unbounded= False;
int numcenters= 0;
int dim= qh hull_dim - 1;
realT dist, offset, angle, zero= 0.0;
midpoint= qh gm_matrix + qh hull_dim * qh hull_dim; /* last row */
for (k= 0; k < dim; k++)
midpoint[k]= (vertex->point[k] + vertexA->point[k])/2;
FOREACHfacet_(centers) {
numcenters++;
if (!facet->visitid)
unbounded= True;
else {
if (!facet->center)
facet->center= qh_facetcenter (facet->vertices);
qh_setappend (&points, facet->center);
}
}
if (numcenters > dim) {
simplex= qh_settemp (qh TEMPsize);
qh_setappend (&simplex, vertex->point);
if (unbounded)
qh_setappend (&simplex, midpoint);
qh_maxsimplex (dim, points, NULL, 0, &simplex);
qh_setdelnth (simplex, 0);
}else if (numcenters == dim) {
if (unbounded)
qh_setappend (&points, midpoint);
simplex= points;
}else {
fprintf(qh ferr, "qh_detvnorm: too few points (%d) to compute separating plane\n", numcenters);
qh_errexit (qh_ERRqhull, NULL, NULL);
}
i= 0;
gmcoord= qh gm_matrix;
point0= SETfirstt_(simplex, pointT);
FOREACHpoint_(simplex) {
if (qh IStracing >= 4)
qh_printmatrix(qh ferr, "qh_detvnorm: Voronoi vertex or midpoint",
&point, 1, dim);
if (point != point0) {
qh gm_row[i++]= gmcoord;
coord= point0;
for (k= dim; k--; )
*(gmcoord++)= *point++ - *coord++;
}
}
qh gm_row[i]= gmcoord; /* does not overlap midpoint, may be used later for qh_areasimplex */
normal= gmcoord;
qh_sethyperplane_gauss (dim, qh gm_row, point0, True,
normal, &offset, &nearzero);
if (qh GOODvertexp == vertexA->point)
inpoint= vertexA->point;
else
inpoint= vertex->point;
zinc_(Zdistio);
dist= qh_distnorm (dim, inpoint, normal, &offset);
if (dist > 0) {
offset= -offset;
normalp= normal;
for (k= dim; k--; ) {
*normalp= -(*normalp);
normalp++;
}
}
if (qh VERIFYoutput || qh PRINTstatistics) {
pointid= qh_pointid (vertex->point);
pointidA= qh_pointid (vertexA->point);
if (!unbounded) {
zinc_(Zdiststat);
dist= qh_distnorm (dim, midpoint, normal, &offset);
if (dist < 0)
dist= -dist;
zzinc_(Zridgemid);
wwmax_(Wridgemidmax, dist);
wwadd_(Wridgemid, dist);
trace4((qh ferr, "qh_detvnorm: points %d %d midpoint dist %2.2g\n",
pointid, pointidA, dist));
for (k= 0; k < dim; k++)
midpoint[k]= vertexA->point[k] - vertex->point[k]; /* overwrites midpoint! */
qh_normalize (midpoint, dim, False);
angle= qh_distnorm (dim, midpoint, normal, &zero); /* qh_detangle uses dim+1 */
if (angle < 0.0)
angle= angle + 1.0;
else
angle= angle - 1.0;
if (angle < 0.0)
angle -= angle;
trace4((qh ferr, "qh_detvnorm: points %d %d angle %2.2g nearzero %d\n",
pointid, pointidA, angle, nearzero));
if (nearzero) {
zzinc_(Zridge0);
wwmax_(Wridge0max, angle);
wwadd_(Wridge0, angle);
}else {
zzinc_(Zridgeok)
wwmax_(Wridgeokmax, angle);
wwadd_(Wridgeok, angle);
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -