📄 geom.c
字号:
/*<html><pre> -<a href="qh-geom.htm"
>-------------------------------</a><a name="TOP">-</a>
geom.c
geometric routines of qhull
see qh-geom.htm and geom.h
copyright (c) 1993-2003 The Geometry Center
infrequent code goes into geom2.c
*/
#include "qhull_a.h"
/*-<a href="qh-geom.htm#TOC"
>-------------------------------</a><a name="distplane">-</a>
qh_distplane( point, facet, dist )
return distance from point to facet
returns:
dist
if qh.RANDOMdist, joggles result
notes:
dist > 0 if point is above facet (i.e., outside)
does not error (for sortfacets)
see:
qh_distnorm in geom2.c
*/
void qh_distplane (pointT *point, facetT *facet, realT *dist) {
coordT *normal= facet->normal, *coordp, randr;
int k;
switch(qh hull_dim){
case 2:
*dist= facet->offset + point[0] * normal[0] + point[1] * normal[1];
break;
case 3:
*dist= facet->offset + point[0] * normal[0] + point[1] * normal[1] + point[2] * normal[2];
break;
case 4:
*dist= facet->offset+point[0]*normal[0]+point[1]*normal[1]+point[2]*normal[2]+point[3]*normal[3];
break;
case 5:
*dist= facet->offset+point[0]*normal[0]+point[1]*normal[1]+point[2]*normal[2]+point[3]*normal[3]+point[4]*normal[4];
break;
case 6:
*dist= facet->offset+point[0]*normal[0]+point[1]*normal[1]+point[2]*normal[2]+point[3]*normal[3]+point[4]*normal[4]+point[5]*normal[5];
break;
case 7:
*dist= facet->offset+point[0]*normal[0]+point[1]*normal[1]+point[2]*normal[2]+point[3]*normal[3]+point[4]*normal[4]+point[5]*normal[5]+point[6]*normal[6];
break;
case 8:
*dist= facet->offset+point[0]*normal[0]+point[1]*normal[1]+point[2]*normal[2]+point[3]*normal[3]+point[4]*normal[4]+point[5]*normal[5]+point[6]*normal[6]+point[7]*normal[7];
break;
default:
*dist= facet->offset;
coordp= point;
for (k= qh hull_dim; k--; )
*dist += *coordp++ * *normal++;
break;
}
zinc_(Zdistplane);
if (!qh RANDOMdist && qh IStracing < 4)
return;
if (qh RANDOMdist) {
randr= qh_RANDOMint;
*dist += (2.0 * randr / qh_RANDOMmax - 1.0) *
qh RANDOMfactor * qh MAXabs_coord;
}
if (qh IStracing >= 4) {
fprintf (qh ferr, "qh_distplane: ");
fprintf (qh ferr, qh_REAL_1, *dist);
fprintf (qh ferr, "from p%d to f%d\n", qh_pointid(point), facet->id);
}
return;
} /* distplane */
/*-<a href="qh-geom.htm#TOC"
>-------------------------------</a><a name="findbest">-</a>
qh_findbest( point, startfacet, bestoutside, qh_ISnewfacets, qh_NOupper, dist, isoutside, numpart )
find facet that is furthest below a point
for upperDelaunay facets
returns facet only if !qh_NOupper and clearly above
input:
starts search at 'startfacet' (can not be flipped)
if !bestoutside (qh_ALL), stops at qh.MINoutside
returns:
best facet (reports error if NULL)
early out if isoutside defined and bestdist > qh.MINoutside
dist is distance to facet
isoutside is true if point is outside of facet
numpart counts the number of distance tests
see also:
qh_findbestnew()
notes:
If merging (testhorizon), searches horizon facets of coplanar best facets because
after qh_distplane, this and qh_partitionpoint are the most expensive in 3-d
avoid calls to distplane, function calls, and real number operations.
caller traces result
Optimized for outside points. Tried recording a search set for qh_findhorizon.
Made code more complicated.
when called by qh_partitionvisible():
indicated by qh_ISnewfacets
qh.newfacet_list is list of simplicial, new facets
qh_findbestnew set if qh_sharpnewfacets returns True (to use qh_findbestnew)
qh.bestfacet_notsharp set if qh_sharpnewfacets returns False
when called by qh_findfacet(), qh_partitionpoint(), qh_partitioncoplanar(),
qh_check_bestdist(), qh_addpoint()
indicated by !qh_ISnewfacets
returns best facet in neighborhood of given facet
this is best facet overall if dist > - qh.MAXcoplanar
or hull has at least a "spherical" curvature
design:
initialize and test for early exit
repeat while there are better facets
for each neighbor of facet
exit if outside facet found
test for better facet
if point is inside and partitioning
test for new facets with a "sharp" intersection
if so, future calls go to qh_findbestnew()
test horizon facets
*/
facetT *qh_findbest (pointT *point, facetT *startfacet,
boolT bestoutside, boolT isnewfacets, boolT noupper,
realT *dist, boolT *isoutside, int *numpart) {
realT bestdist= -REALmax/2 /* avoid underflow */;
facetT *facet, *neighbor, **neighborp;
facetT *bestfacet= NULL, *lastfacet= NULL;
int oldtrace= qh IStracing;
unsigned int visitid= ++qh visit_id;
int numpartnew=0;
boolT testhorizon = True; /* needed if precise, e.g., rbox c D6 | qhull Q0 Tv */
zinc_(Zfindbest);
if (qh IStracing >= 3 || (qh TRACElevel && qh TRACEpoint >= 0 && qh TRACEpoint == qh_pointid (point))) {
if (qh TRACElevel > qh IStracing)
qh IStracing= qh TRACElevel;
fprintf (qh ferr, "qh_findbest: point p%d starting at f%d isnewfacets? %d, unless %d exit if > %2.2g\n",
qh_pointid(point), startfacet->id, isnewfacets, bestoutside, qh MINoutside);
fprintf(qh ferr, " testhorizon? %d noupper? %d", testhorizon, noupper);
fprintf (qh ferr, " Last point added was p%d.", qh furthest_id);
fprintf(qh ferr, " Last merge was #%d. max_outside %2.2g\n", zzval_(Ztotmerge), qh max_outside);
}
if (isoutside)
*isoutside= True;
if (!startfacet->flipped) { /* test startfacet */
*numpart= 1;
qh_distplane (point, startfacet, dist); /* this code is duplicated below */
if (!bestoutside && *dist >= qh MINoutside
&& (!startfacet->upperdelaunay || !noupper)) {
bestfacet= startfacet;
goto LABELreturn_best;
}
bestdist= *dist;
if (!startfacet->upperdelaunay) {
bestfacet= startfacet;
}
}else
*numpart= 0;
startfacet->visitid= visitid;
facet= startfacet;
while (facet) {
trace4((qh ferr, "qh_findbest: neighbors of f%d, bestdist %2.2g f%d\n",
facet->id, bestdist, getid_(bestfacet)));
lastfacet= facet;
FOREACHneighbor_(facet) {
if (!neighbor->newfacet && isnewfacets)
continue;
if (neighbor->visitid == visitid)
continue;
neighbor->visitid= visitid;
if (!neighbor->flipped) { /* code duplicated above */
(*numpart)++;
qh_distplane (point, neighbor, dist);
if (*dist > bestdist) {
if (!bestoutside && *dist >= qh MINoutside
&& (!neighbor->upperdelaunay || !noupper)) {
bestfacet= neighbor;
goto LABELreturn_best;
}
if (!neighbor->upperdelaunay) {
bestfacet= neighbor;
bestdist= *dist;
break; /* switch to neighbor */
}else if (!bestfacet) {
bestdist= *dist;
break; /* switch to neighbor */
}
} /* end of *dist>bestdist */
} /* end of !flipped */
} /* end of FOREACHneighbor */
facet= neighbor; /* non-NULL only if *dist>bestdist */
} /* end of while facet (directed search) */
if (isnewfacets) {
if (!bestfacet) {
bestdist= -REALmax/2;
bestfacet= qh_findbestnew (point, startfacet->next, &bestdist, bestoutside, isoutside, &numpartnew);
testhorizon= False; /* qh_findbestnew calls qh_findbesthorizon */
}else if (!qh findbest_notsharp && bestdist < - qh DISTround) {
if (qh_sharpnewfacets()) {
/* seldom used, qh_findbestnew will retest all facets */
zinc_(Zfindnewsharp);
bestfacet= qh_findbestnew (point, bestfacet, &bestdist, bestoutside, isoutside, &numpartnew);
testhorizon= False; /* qh_findbestnew calls qh_findbesthorizon */
qh findbestnew= True;
}else
qh findbest_notsharp= True;
}
}
if (!bestfacet)
bestfacet= qh_findbestlower (lastfacet, point, &bestdist, numpart);
if (testhorizon)
bestfacet= qh_findbesthorizon (!qh_IScheckmax, point, bestfacet, noupper, &bestdist, &numpartnew);
*dist= bestdist;
if (isoutside && bestdist < qh MINoutside)
*isoutside= False;
LABELreturn_best:
zadd_(Zfindbesttot, *numpart);
zmax_(Zfindbestmax, *numpart);
(*numpart) += numpartnew;
qh IStracing= oldtrace;
return bestfacet;
} /* findbest */
/*-<a href="qh-geom.htm#TOC"
>-------------------------------</a><a name="findbesthorizon">-</a>
qh_findbesthorizon( qh_IScheckmax, point, startfacet, qh_NOupper, &bestdist, &numpart )
search coplanar and better horizon facets from startfacet/bestdist
ischeckmax turns off statistics and minsearch update
all arguments must be initialized
returns (ischeckmax):
best facet
returns (!ischeckmax):
best facet that is not upperdelaunay
allows upperdelaunay that is clearly outside
returns:
bestdist is distance to bestfacet
numpart -- updates number of distance tests
notes:
no early out -- use qh_findbest() or qh_findbestnew()
Searches coplanar or better horizon facets
when called by qh_check_maxout() (qh_IScheckmax)
startfacet must be closest to the point
Otherwise, if point is beyond and below startfacet, startfacet may be a local minimum
even though other facets are below the point.
updates facet->maxoutside for good, visited facets
may return NULL
searchdist is qh.max_outside + 2 * DISTround
+ max( MINvisible('Vn'), MAXcoplanar('Un'));
This setting is a guess. It must be at least max_outside + 2*DISTround
because a facet may have a geometric neighbor across a vertex
design:
for each horizon facet of coplanar best facets
continue if clearly inside
unless upperdelaunay or clearly outside
update best facet
*/
facetT *qh_findbesthorizon (boolT ischeckmax, pointT* point, facetT *startfacet, boolT noupper, realT *bestdist, int *numpart) {
facetT *bestfacet= startfacet;
realT dist;
facetT *neighbor, **neighborp, *facet;
facetT *nextfacet= NULL; /* optimize last facet of coplanarset */
int numpartinit= *numpart, coplanarset_size;
unsigned int visitid= ++qh visit_id;
boolT newbest= False; /* for tracing */
realT minsearch, searchdist; /* skip facets that are too far from point */
if (!ischeckmax) {
zinc_(Zfindhorizon);
}else {
#if qh_MAXoutside
if ((!qh ONLYgood || startfacet->good) && *bestdist > startfacet->maxoutside)
startfacet->maxoutside= *bestdist;
#endif
}
searchdist= qh_SEARCHdist; /* multiple of qh.max_outside and precision constants */
minsearch= *bestdist - searchdist;
if (ischeckmax) {
/* Always check coplanar facets. Needed for RBOX 1000 s Z1 G1e-13 t996564279 | QHULL Tv */
minimize_(minsearch, -searchdist);
}
coplanarset_size= 0;
facet= startfacet;
while (True) {
trace4((qh ferr, "qh_findbesthorizon: neighbors of f%d bestdist %2.2g f%d ischeckmax? %d noupper? %d minsearch %2.2g searchdist %2.2g\n",
facet->id, *bestdist, getid_(bestfacet), ischeckmax, noupper,
minsearch, searchdist));
FOREACHneighbor_(facet) {
if (neighbor->visitid == visitid)
continue;
neighbor->visitid= visitid;
if (!neighbor->flipped) {
qh_distplane (point, neighbor, &dist);
(*numpart)++;
if (dist > *bestdist) {
if (!neighbor->upperdelaunay || ischeckmax || (!noupper && dist >= qh MINoutside)) {
bestfacet= neighbor;
*bestdist= dist;
newbest= True;
if (!ischeckmax) {
minsearch= dist - searchdist;
if (dist > *bestdist + searchdist) {
zinc_(Zfindjump); /* everything in qh.coplanarset at least searchdist below */
coplanarset_size= 0;
}
}
}
}else if (dist < minsearch)
continue; /* if ischeckmax, dist can't be positive */
#if qh_MAXoutside
if (ischeckmax && dist > neighbor->maxoutside)
neighbor->maxoutside= dist;
#endif
} /* end of !flipped */
if (nextfacet) {
if (!coplanarset_size++) {
SETfirst_(qh coplanarset)= nextfacet;
SETtruncate_(qh coplanarset, 1);
}else
qh_setappend (&qh coplanarset, nextfacet); /* Was needed for RBOX 1000 s W1e-13 P0 t996547055 | QHULL d Qbb Qc Tv
and RBOX 1000 s Z1 G1e-13 t996564279 | qhull Tv */
}
nextfacet= neighbor;
} /* end of EACHneighbor */
facet= nextfacet;
if (facet)
nextfacet= NULL;
else if (!coplanarset_size)
break;
else if (!--coplanarset_size) {
facet= SETfirst_(qh coplanarset);
SETtruncate_(qh coplanarset, 0);
}else
facet= (facetT*)qh_setdellast (qh coplanarset);
} /* while True, for each facet in qh.coplanarset */
if (!ischeckmax) {
zadd_(Zfindhorizontot, *numpart - numpartinit);
zmax_(Zfindhorizonmax, *numpart - numpartinit);
if (newbest)
zinc_(Zparthorizon);
}
trace4((qh ferr, "qh_findbesthorizon: newbest? %d bestfacet f%d bestdist %2.2g\n", newbest, getid_(bestfacet), *bestdist));
return bestfacet;
} /* findbesthorizon */
/*-<a href="qh-geom.htm#TOC"
>-------------------------------</a><a name="findbestnew">-</a>
qh_findbestnew( point, startfacet, dist, isoutside, numpart )
find best newfacet for point
searches all of qh.newfacet_list starting at startfacet
searches horizon facets of coplanar best newfacets
searches all facets if startfacet == qh.facet_list
returns:
best new or horizon facet that is not upperdelaunay
early out if isoutside and not 'Qf'
dist is distance to facet
isoutside is true if point is outside of facet
numpart is number of distance tests
notes:
Always used for merged new facets (see qh_USEfindbestnew)
Avoids upperdelaunay facet unless (isoutside and outside)
Uses qh.visit_id, qh.coplanarset.
If share visit_id with qh_findbest, coplanarset is incorrect.
If merging (testhorizon), searches horizon facets of coplanar best facets because
a point maybe coplanar to the bestfacet, below its horizon facet,
and above a horizon facet of a coplanar newfacet. For example,
rbox 1000 s Z1 G1e-13 | qhull
rbox 1000 s W1e-13 P0 t992110337 | QHULL d Qbb Qc
qh_findbestnew() used if
qh_sharpnewfacets -- newfacets contains a sharp angle
if many merges, qh_premerge found a merge, or 'Qf' (qh.findbestnew)
see also:
qh_partitionall() and qh_findbest()
design:
for each new facet starting from startfacet
test distance from point to facet
return facet if clearly outside
unless upperdelaunay and a lowerdelaunay exists
update best facet
test horizon facets
*/
facetT *qh_findbestnew (pointT *point, facetT *startfacet,
realT *dist, boolT bestoutside, boolT *isoutside, int *numpart) {
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -