📄 qhull.c
字号:
/*<html><pre> -<a href="qh-c.htm#qhull"
>-------------------------------</a><a name="TOP">-</a>
qhull.c
Quickhull algorithm for convex hulls
qhull() and top-level routines
see qh-c.htm, qhull.h, unix.c
see qhull_a.h for internal functions
copyright (c) 1993-1999 The Geometry Center
*/
#include "qhull_a.h"
/*============= functions in alphabetic order after qhull() =======*/
/*-<a href="qh-c.htm#qhull"
>-------------------------------</a><a name="qhull">-</a>
qh_qhull()
compute DIM3 convex hull of qh.num_points starting at qh.first_point
qh contains all global options and variables
returns:
returns polyhedron
qh.facet_list, qh.num_facets, qh.vertex_list, qh.num_vertices,
returns global variables
qh.hulltime, qh.max_outside, qh.interior_point, qh.max_vertex, qh.min_vertex
returns precision constants
qh.ANGLEround, centrum_radius, cos_max, DISTround, MAXabs_coord, ONEmerge
notes:
unless needed for output
qh.max_vertex and qh.min_vertex are max/min due to merges
see:
to add individual points to either qh.num_points
use qh_addpoint()
if qh.GETarea
qh_produceoutput() returns qh.totarea and qh.totvol via qh_getarea()
design:
record starting time
initialize hull and partition points
build convex hull
unless early termination
update facet->maxoutside for vertices, coplanar, and near-inside points
error if temporary sets exist
record end time
*/
void qh_qhull (void) {
int numoutside;
qh hulltime= qh_CPUclock;
if (qh RERUN || qh JOGGLEmax < REALmax/2)
qh_build_withrestart();
else {
qh_initbuild();
qh_buildhull();
}
if (!qh STOPpoint && !qh STOPcone) {
if (qh ZEROall_ok && !qh TESTvneighbors && qh MERGEexact)
qh_checkzero( qh_ALL);
if (qh ZEROall_ok && !qh TESTvneighbors && !qh WAScoplanar) {
trace2((qh ferr, "qh_qhull: all facets are clearly convex and no coplanar points. Post-merging and check of maxout not needed.\n"));
}else {
if (qh MERGEexact || (qh hull_dim > qh_DIMreduceBuild && qh PREmerge))
qh_postmerge ("First post-merge", qh premerge_centrum, qh premerge_cos,
(qh POSTmerge ? False : qh TESTvneighbors));
else if (!qh POSTmerge && qh TESTvneighbors)
qh_postmerge ("For testing vertex neighbors", qh premerge_centrum,
qh premerge_cos, True);
if (qh POSTmerge)
qh_postmerge ("For post-merging", qh postmerge_centrum,
qh postmerge_cos, qh TESTvneighbors);
if (qh visible_list == qh facet_list) { /* i.e., merging done */
qh findbestnew= True;
qh_partitionvisible (/*visible_list, newfacet_list*/ !qh_ALL, &numoutside);
qh findbestnew= False;
qh_deletevisible (/*qh visible_list*/);
qh_resetlists (False /*qh visible_list newvertex_list newfacet_list */);
}
if (qh DOcheckmax){
if (qh REPORTfreq) {
qh_buildtracing (NULL, NULL);
fprintf (qh ferr, "\nTesting all coplanar points.\n");
}
qh_check_maxout();
}
}
}
if (qh KEEPnearinside && !qh maxoutdone)
qh_nearcoplanar();
if (qh_setsize ((setT*)qhmem.tempstack) != 0) {
fprintf (qh ferr, "qhull internal error (qh_qhull): temporary sets not empty (%d)\n",
qh_setsize ((setT*)qhmem.tempstack));
qh_errexit (qh_ERRqhull, NULL, NULL);
}
qh hulltime= qh_CPUclock - qh hulltime;
qh QHULLfinished= True;
trace1((qh ferr, "qh_qhull: algorithm completed\n"));
} /* qhull */
/*-<a href="qh-c.htm#qhull"
>-------------------------------</a><a name="addpoint">-</a>
qh_addpoint( furthest, facet, checkdist )
add point (usually furthest point) above facet to hull
if checkdist,
check that point is above facet.
if point is not outside of the hull, uses qh_partitioncoplanar()
else if facet specified,
assumes that point is above facet (major damage if below)
for Delaunay triangulations,
Use qh_setdelaunay() to lift point to paraboloid and scale by 'Qbb' if needed
Do not use options 'Qbk', 'QBk', or 'QbB' since they scale the coordinates.
returns:
returns False if user requested an early termination
qh.visible_list, newfacet_list, delvertex_list, NEWfacets may be defined
updates qh.facet_list, qh.num_facets, qh.vertex_list, qh.num_vertices
clear qh.maxoutdone (will need to call qh_check_maxout() for facet->maxoutside)
if unknown point, adds a pointer to qh.other_points
do not deallocate the point's coordinates
notes:
uses qh.visible_list, qh.newfacet_list, qh.delvertex_list, qh.NEWfacets
design:
check point in qh.first_point/.num_points
if checkdist
if point not above facet
partition coplanar point
exit
exit if pre STOPpoint requested
find horizon and visible facets for point
make new facets for point to horizon
make hyperplanes for point
compute balance statistics
match neighboring new facets
update vertex neighbors and delete interior vertices
exit if STOPcone requested
merge non-convex new facets
check for using qh_findbestnew() instead of qh_findbest()
partition outside points from visible facets
delete visible facets
check polyhedron if requested
exit if post STOPpoint requested
reset working lists of facets and vertices
*/
boolT qh_addpoint (pointT *furthest, facetT *facet, boolT checkdist) {
int goodvisible, goodhorizon;
vertexT *vertex;
facetT *newfacet;
realT dist, newbalance, pbalance;
boolT isoutside= False;
int numpart, numpoints, numnew, firstnew;
qh maxoutdone= False;
if (qh_pointid (furthest) == -1)
qh_setappend (&qh other_points, furthest);
if (!facet) {
fprintf (qh ferr, "qh_addpoint: NULL facet. Use qh_findbestfacet\n");
qh_errexit (qh_ERRqhull, NULL, NULL);
}
if (checkdist) {
facet= qh_findbest (furthest, facet, !qh_ALL, False, !qh_NOupper,
&dist, &isoutside, &numpart);
zzadd_(Zpartition, numpart);
if (!isoutside) {
zinc_(Znotmax); /* last point of outsideset is no longer furthest. */
facet->notfurthest= True;
qh_partitioncoplanar (furthest, facet, &dist);
return True;
}
}
qh_buildtracing (furthest, facet);
if (qh STOPpoint < 0 && qh furthest_id == -qh STOPpoint-1) {
facet->notfurthest= True;
return False;
}
qh_findhorizon (furthest, facet, &goodvisible, &goodhorizon);
if (qh ONLYgood && !(goodvisible+goodhorizon) && !qh GOODclosest) {
zinc_(Znotgood);
facet->notfurthest= True;
/* last point of outsideset is no longer furthest. This is ok
since all points of the outside are likely to be bad */
qh_resetlists (False /*qh visible_list newvertex_list newfacet_list */);
return True;
}
zzinc_(Zprocessed);
firstnew= qh facet_id;
vertex= qh_makenewfacets (furthest /*visible_list, attaches if !ONLYgood */);
qh_makenewplanes (/* newfacet_list */);
numnew= qh facet_id - firstnew;
newbalance= numnew - (realT) (qh num_facets-qh num_visible)
* qh hull_dim/qh num_vertices;
wadd_(Wnewbalance, newbalance);
wadd_(Wnewbalance2, newbalance * newbalance);
if (qh ONLYgood
&& !qh_findgood (qh newfacet_list, goodhorizon) && !qh GOODclosest) {
FORALLnew_facets
qh_delfacet (newfacet);
qh_delvertex (vertex);
qh_resetlists (True /*qh visible_list newvertex_list newfacet_list */);
zinc_(Znotgoodnew);
facet->notfurthest= True;
return True;
}
if (qh ONLYgood)
qh_attachnewfacets(/*visible_list*/);
qh_matchnewfacets();
qh_updatevertices();
if (qh STOPcone && qh furthest_id == qh STOPcone-1) {
facet->notfurthest= True;
return False; /* visible_list etc. still defined */
}
if (qh PREmerge || qh MERGEexact) {
qh_premerge (vertex, qh premerge_centrum, qh premerge_cos);
if (zzval_(Ztotmerge) > qh_USEfindbestnew)
qh findbestnew= True;
else {
FORALLnew_facets {
if (!newfacet->simplicial) {
qh findbestnew= True; /* use qh_findbestnew instead of qh_findbest*/
break;
}
}
}
}else if (qh BESToutside)
qh findbestnew= True;
qh_partitionvisible (/*visible_list, newfacet_list*/ !qh_ALL, &numpoints);
qh findbestnew= False;
qh findbest_notsharp= False;
zinc_(Zpbalance);
pbalance= numpoints - (realT) qh hull_dim /* assumes all points extreme */
* (qh num_points - qh num_vertices)/qh num_vertices;
wadd_(Wpbalance, pbalance);
wadd_(Wpbalance2, pbalance * pbalance);
qh_deletevisible (/*qh visible_list*/);
zmax_(Zmaxvertex, qh num_vertices);
qh NEWfacets= False;
if (qh IStracing >= 4)
qh_printfacetlist (qh newfacet_list, NULL, True);
if (qh CHECKfrequently) {
if (qh num_facets < 50)
qh_checkpolygon (qh facet_list);
else
qh_checkpolygon (qh newfacet_list);
}
if (qh STOPpoint > 0 && qh furthest_id == qh STOPpoint-1)
return False;
qh_resetlists (True /*qh visible_list newvertex_list newfacet_list */);
trace2((qh ferr, "qh_addpoint: added p%d new facets %d new balance %2.2g point balance %2.2g\n",
qh_pointid (furthest), numnew, newbalance, pbalance));
if (qh hull_dim > 3 && qh TRACEpoint == qh_pointid (furthest)) {
qh IStracing= 0;
qhmem.IStracing= 0;
}
return True;
} /* addpoint */
/*-<a href="qh-c.htm#qhull"
>-------------------------------</a><a name="build_withrestart">-</a>
qh_build_withrestart()
allow restarts due to qh.JOGGLEmax while calling qh_buildhull()
qh.FIRSTpoint/qh.NUMpoints is point array
it may be moved by qh_joggleinput()
*/
void qh_build_withrestart (void) {
int restart;
qh ALLOWrestart= True;
while (True) {
restart= setjmp (qh restartexit); /* simple statement for CRAY J916 */
if (restart) { /* only from qh_precision() */
zzinc_(Zretry);
wmax_(Wretrymax, qh JOGGLEmax);
qh ERREXITcalled= False;
qh STOPcone= True; /* if break, prevents normal output */
}
if (!qh RERUN && qh JOGGLEmax < REALmax/2) {
if (qh build_cnt > qh_JOGGLEmaxretry) {
fprintf(qh ferr, "\n\
qhull precision error: %d attempts to construct a convex hull\n\
with joggled input. Increase joggle above 'QJ%2.2g'\n\
or modify qh_JOGGLE... parameters in user.h\n",
qh build_cnt, qh JOGGLEmax);
qh_errexit (qh_ERRqhull, NULL, NULL);
}
if (qh build_cnt && !restart)
break;
}else if (qh build_cnt && qh build_cnt >= qh RERUN)
break;
qh STOPcone= False;
qh_freebuild (True); /* first call is a nop */
qh build_cnt++;
if (!qh qhull_optionsiz)
qh qhull_optionsiz= strlen (qh qhull_options);
else {
qh qhull_options [qh qhull_optionsiz]= '\0';
qh qhull_optionlen= 80;
}
qh_option("_run", &qh build_cnt, NULL);
if (qh build_cnt == qh RERUN) {
qh IStracing= qh TRACElastrun; /* duplicated from qh_initqhull_globals */
if (qh TRACEpoint != -1 || qh TRACEdist < REALmax/2 || qh TRACEmerge) {
qh TRACElevel= (qh IStracing? qh IStracing : 3);
qh IStracing= 0;
}
qhmem.IStracing= qh IStracing;
}
if (qh JOGGLEmax < REALmax/2)
qh_joggleinput();
qh_initbuild();
qh_buildhull();
if (qh JOGGLEmax < REALmax/2 && !qh MERGING)
qh_checkconvex (qh facet_list, qh_ALGORITHMfault);
}
qh ALLOWrestart= False;
} /* qh_build_withrestart */
/*-<a href="qh-c.htm#qhull"
>-------------------------------</a><a name="buildhull">-</a>
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -