📄 geom.c
字号:
/*<html><pre> -<a href="qh-c.htm"
>-------------------------------</a><a name="TOP">-</a>
geom.c
geometric routines of qhull
see qh-c.htm and geom.h
copyright (c) 1993-1999 The Geometry Center
infrequent code goes into geom2.c
*/
#include "qhull_a.h"
/*-<a href="qh-c.htm#geom"
>-------------------------------</a><a name="backnormal">-</a>
qh_backnormal( rows, numrow, numcol, sign, normal, nearzero )
given an upper-triangular rows array and a sign,
solve for normal equation x using back substitution over rows U
returns:
normal= x
if will not be able to divzero() when normalized (qh.MINdenom_2 and qh.MINdenom_1_2),
if fails on last row
this means that the hyperplane intersects [0,..,1]
sets last coordinate of normal to sign
otherwise
sets tail of normal to [...,sign,0,...], i.e., solves for b= [0...0]
sets nearzero
notes:
assumes numrow == numcol-1
see Golub & van Loan 4.4-9 for back substitution
solves Ux=b where Ax=b and PA=LU
b= [0,...,0,sign or 0] (sign is either -1 or +1)
last row of A= [0,...,0,1]
1) Ly=Pb == y=b since P only permutes the 0's of b
design:
for each row from end
perform back substitution
if near zero
use qh_divzero for division
if zero divide and not last row
set tail of normal to 0
*/
void qh_backnormal (realT **rows, int numrow, int numcol, boolT sign,
coordT *normal, boolT *nearzero) {
int i, j;
coordT *normalp, *normal_tail, *ai, *ak;
realT diagonal;
boolT waszero;
int zerocol= -1;
normalp= normal + numcol - 1;
*normalp--= (sign ? -1.0 : 1.0);
for(i= numrow; i--; ) {
*normalp= 0.0;
ai= rows[i] + i + 1;
ak= normalp+1;
for(j= i+1; j < numcol; j++)
*normalp -= *ai++ * *ak++;
diagonal= (rows[i])[i];
if (fabs_(diagonal) > qh MINdenom_2)
*(normalp--) /= diagonal;
else {
waszero= False;
*normalp= qh_divzero (*normalp, diagonal, qh MINdenom_1_2, &waszero);
if (waszero) {
zerocol= i;
*(normalp--)= (sign ? -1.0 : 1.0);
for (normal_tail= normalp+2; normal_tail < normal + numcol; normal_tail++)
*normal_tail= 0.0;
}else
normalp--;
}
}
if (zerocol != -1) {
zzinc_(Zback0);
*nearzero= True;
trace4((qh ferr, "qh_backnormal: zero diagonal at column %d.\n", i));
qh_precision ("zero diagonal on back substitution");
}
} /* backnormal */
/*-<a href="qh-c.htm#geom"
>-------------------------------</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-c.htm#geom"
>-------------------------------</a><a name="findbest">-</a>
qh_findbest( point, startfacet, bestoutside, newfacts, dist, isoutside, numpart )
find facet that is furthest below a point
searches neighbors of coplanar and most flipped facets
for upperDelaunay facets
searches if coplanar (within searchdist)
returns facet only if !noupper and clearly above
input:
starts search at 'startfacet' (can not be flipped)
if !bestoutside, stops at qh.MINoutside (DISTround if precise)
returns:
best facet
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:
uses qh.visit_id, qh.searchset
caller traces the results
after qh_distplane, this and qh_partitionpoint are the most expensive in 3-d
avoid calls to distplane, function calls and real number operations.
when called by qh_partitionvisible():
indicated by newfacets and isoutside defined
qh.newfacet_list is list of simplicial, new facets
qh_findbestnew set if qh_findbestsharp returns True (to use qh_findbestnew)
qh.bestfacet_notsharp set if qh_findbestsharp returns False
searches horizon of best facet unless "exact" and !bestoutside
searchdist is 2 * DISTround
when called by qh_check_maxout()
indicated by bestoutside and !newfacets and isoutside == NULL
startfacet must be closest to 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
when called by findfacet() and check_bestdist()
indicated by !newfacets and isoutside defined
searchdist is same as qh_check_maxout()
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 facets to search (searchset)
for each neighbor of facet
exit if outside facet found
restart searchset if neighbor is much better
append flipped and nearby facets to searchset
if point is inside and partitioning
test for new facets with a "sharp" intersection
if so, future calls go to qh_findbestsharp
if testhorizon
also test facet's horizon facet
*/
facetT *qh_findbest (pointT *point, facetT *startfacet,
boolT bestoutside, boolT newfacets, boolT noupper,
realT *dist, boolT *isoutside, int *numpart) {
realT bestdist= -REALmax/2 /* avoid underflow */, searchdist;
realT cutoff, mincutoff; /* skip facets that are too far from point */
facetT *facet, *neighbor, **neighborp, *bestfacet= NULL;
int oldtrace= qh IStracing;
int searchsize= 0; /* non-zero if searchset defined */
boolT newbest;
boolT ischeckmax= bestoutside && !newfacets && !isoutside;
boolT ispartition= newfacets && isoutside;
boolT isfindfacet= !newfacets && isoutside;
boolT testhorizon = ispartition && (bestoutside || qh APPROXhull || qh MERGING);
if (!ischeckmax && !ispartition && !isfindfacet) {
fprintf (qh ferr, "qhull internal error (qh_findbest): unknown combination of arguments\n");
qh_errexit (qh_ERRqhull, startfacet, NULL);
}
if (qh TRACElevel && qh TRACEpoint >= 0 && qh TRACEpoint == qh_pointid (point)) {
qh IStracing= qh TRACElevel;
fprintf (qh ferr, "qh_findbest: point p%d starting at f%d bestoutside? %d newfacets %d\n",
qh TRACEpoint, startfacet->id, bestoutside, newfacets);
fprintf(qh ferr, " ischeckmax %d ispartition %d isfindfacet %d testhorizon %d\n",
ischeckmax, ispartition, isfindfacet, testhorizon);
fprintf (qh ferr, " Last point added to hull was p%d.", qh furthest_id);
fprintf(qh ferr, " Last merge was #%d.\n", zzval_(Ztotmerge));
}
if (isoutside)
*isoutside= True;
if (!startfacet->flipped) {
*numpart= 1;
qh_distplane (point, startfacet, dist); /* this code is duplicated below */
if (!startfacet->upperdelaunay || (!noupper && *dist >= qh MINoutside)) {
bestdist= *dist;
bestfacet= startfacet;
if (!bestoutside && *dist >= qh MINoutside)
goto LABELreturn_best;
}
#if qh_MAXoutside
if (ischeckmax && (!qh ONLYgood || startfacet->good) && *dist > startfacet->maxoutside)
startfacet->maxoutside= *dist;
#endif
}
if (ispartition)
searchdist= 2 * qh DISTround;
else
searchdist= qh max_outside + 2 * qh DISTround
+ fmax_( qh MINvisible, qh MAXcoplanar);
cutoff= bestdist - searchdist;
mincutoff= 0;
if (ischeckmax) {
mincutoff= -(qh DISTround - fmax_(qh MINvisible, qh MAXcoplanar));
minimize_(cutoff, mincutoff);
}
startfacet->visitid= ++qh visit_id;
facet= startfacet;
do { /* search neighbors of coplanar, upperdelaunay, and flipped facets */
if (True) {
LABELrestart: /* directed search whenever improvement > searchdist */
newbest= False;
trace4((qh ferr, "qh_findbest: neighbors of f%d, bestdist %2.2g cutoff %2.2g searchdist %2.2g\n",
facet->id, bestdist, cutoff, searchdist));
FOREACHneighbor_(facet) {
if (ispartition && !neighbor->newfacet)
continue;
if (!neighbor->flipped) {
if (neighbor->visitid == qh visit_id)
continue;
neighbor->visitid= qh visit_id;
(*numpart)++;
qh_distplane (point, neighbor, dist);
if (!bestoutside && *dist >= qh MINoutside
&& (!noupper || !facet->upperdelaunay)) {
bestfacet= neighbor;
goto LABELreturn_best;
}
#if qh_MAXoutside
if (ischeckmax) {
if ((!qh ONLYgood || neighbor->good)
&& *dist > neighbor->maxoutside)
neighbor->maxoutside= *dist;
else if (bestfacet && *dist < cutoff)
continue;
}else
#endif /* dangling else! */
if (bestfacet && *dist < cutoff)
continue;
if (*dist > bestdist) {
if (!neighbor->upperdelaunay
|| (bestoutside && !noupper && *dist >= qh MINoutside)) {
if (ischeckmax && qh_MAXoutside) {
bestdist= *dist;
bestfacet= neighbor;
cutoff= bestdist - searchdist;
minimize_(cutoff, mincutoff);
}else if (*dist > bestdist + searchdist) {
bestdist= *dist;
bestfacet= neighbor;
cutoff= bestdist - searchdist;
searchsize= 0;
facet= neighbor;
if (newbest) /* newbest may be coplanar with facet */
facet->visitid= ++qh visit_id;
goto LABELrestart; /* repeat with a new facet */
}else {
bestdist= *dist;
bestfacet= neighbor;
cutoff= bestdist - searchdist;
}
newbest= True;
}
}
}
if (!searchsize++) {
SETfirst_(qh searchset) = neighbor;
qh_settruncate (qh searchset, 1);
}else
qh_setappend (&qh searchset, neighbor);
} /* FOREACHneighbor */
} /* while (restart) */
}while
(searchsize && (facet= (facetT*)qh_setdellast (qh searchset)));
if (!ischeckmax) {
if (!bestfacet) {
fprintf (qh ferr, "qh_findbest: point p%d starting at f%d bestoutside? %d newfacets %d\n",
qh TRACEpoint, startfacet->id, bestoutside, newfacets);
fprintf(qh ferr, "\n\
qh_findbest: all neighbors of facet %d are flipped or upper Delaunay.\n\
Please report this error to qhull_bug@geom.umn.edu with the input and all of the output.\n",
startfacet->id);
qh FORCEoutput= True;
qh_errexit (qh_ERRqhull, startfacet, NULL);
}
if (ispartition && !qh findbest_notsharp && bestdist < - qh DISTround) {
if (qh_findbestsharp ( point, &bestfacet, &bestdist, numpart))
qh findbestnew= True;
else
qh findbest_notsharp= True;
}
if (testhorizon) {
facet= SETfirstt_(bestfacet->neighbors, facetT);
trace4((qh ferr, "qh_findbest: horizon facet f%d\n", facet->id));
(*numpart)++;
qh_distplane (point, facet, dist);
if (*dist > bestdist
&& (!facet->upperdelaunay || (!noupper && *dist >= qh MINoutside))) {
bestdist= *dist;
bestfacet= facet;
}
}
}
*dist= bestdist;
if (isoutside && bestdist < qh MINoutside)
*isoutside= False;
LABELreturn_best:
qh IStracing= oldtrace;
return bestfacet;
} /* findbest */
/*-<a href="qh-c.htm#geom"
>-------------------------------</a><a name="findbestnew">-</a>
qh_findbestnew( point, startfacet, dist, isoutside, numpart )
find best newfacet for point
searches new facets starting at startfacet
returns:
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -