📄 qhull_geom2.cxx
字号:
/*<html><pre> -<a href="qh-c.htm#geom"
>-------------------------------</a><a name="TOP">-</a>
geom2.c
infrequently used geometric routines of qhull
see qh-c.htm and geom.h
copyright (c) 1993-1999 The Geometry Center
frequently used code goes into geom.c
*/
#include <ivp_physics.hxx>
#include "qhull_a.hxx"
/*================== functions in alphabetic order ============*/
/*-<a href="qh-c.htm#geom"
>-------------------------------</a><a name="copypoints">-</a>
qh_copypoints( points, numpoints, dimension)
return malloc'd copy of points
*/
coordT *qh_copypoints (coordT *points, int numpoints, int dimension) {
int size;
coordT *newpoints;
size= numpoints * dimension * sizeof(coordT);
if (!(newpoints=(coordT*)malloc(size))) {
fprintf(qh ferr, "qhull error: insufficient memory to copy %d points\n",
numpoints);
qh_errexit(qh_ERRmem, NULL, NULL);
}
memcpy ((char *)newpoints, (char *)points, size);
return newpoints;
} /* copypoints */
/*-<a href="qh-c.htm#geom"
>-------------------------------</a><a name="crossproduct">-</a>
qh_crossproduct( dim, vecA, vecB, vecC )
crossproduct of 2 dim vectors
C= A x B
notes:
from Glasner, Graphics Gems I, p. 639
only defined for dim==3
*/
void qh_crossproduct (int dim, realT vecA[3], realT vecB[3], realT vecC[3]){
if (dim == 3) {
vecC[0]= det2_(vecA[1], vecA[2],
vecB[1], vecB[2]);
vecC[1]= - det2_(vecA[0], vecA[2],
vecB[0], vecB[2]);
vecC[2]= det2_(vecA[0], vecA[1],
vecB[0], vecB[1]);
}
} /* vcross */
/*-<a href="qh-c.htm#geom"
>-------------------------------</a><a name="determinant">-</a>
qh_determinant( rows, dim, nearzero )
compute signed determinant of a square matrix
uses qh.NEARzero to test for degenerate matrices
returns:
determinant
overwrites rows and the matrix
if dim == 2 or 3
nearzero iff determinant < qh NEARzero[dim-1]
(not quite correct, not critical)
if dim >= 4
nearzero iff diagonal[k] < qh NEARzero[k]
*/
realT qh_determinant (realT **rows, int dim, boolT *nearzero) {
realT det=0;
int i;
boolT sign= False;
*nearzero= False;
if (dim < 2) {
ivp_message( "qhull internal error (qh_determinate): only implemented for dimension >= 2\n");
qh_errexit (qh_ERRqhull, NULL, NULL);
}else if (dim == 2) {
det= det2_(rows[0][0], rows[0][1],
rows[1][0], rows[1][1]);
if (fabs_(det) < qh NEARzero[1]) /* not really correct, what should this be? */
*nearzero= True;
}else if (dim == 3) {
det= det3_(rows[0][0], rows[0][1], rows[0][2],
rows[1][0], rows[1][1], rows[1][2],
rows[2][0], rows[2][1], rows[2][2]);
if (fabs_(det) < qh NEARzero[2]) /* not really correct, what should this be? */
*nearzero= True;
}else {
qh_gausselim(rows, dim, dim, &sign, nearzero); /* if nearzero, diagonal still ok*/
det= 1.0;
for (i= dim; i--; )
det *= (rows[i])[i];
if (sign)
det= -det;
}
return det;
} /* determinant */
/*-<a href="qh-c.htm#geom"
>-------------------------------</a><a name="detjoggle">-</a>
qh_detjoggle( points, numpoints, dimension )
determine default max joggle for point array
as qh_distround * qh_JOGGLEdefault
returns:
initial value for JOGGLEmax from points and REALepsilon
notes:
computes DISTround since qh_maxmin not called yet
if qh SCALElast, last dimension will be scaled later to MAXwidth
loop duplicated from qh_maxmin
*/
realT qh_detjoggle (pointT *points, int numpoints, int dimension) {
realT abscoord, distround, joggle, maxcoord, mincoord;
pointT *point, *pointtemp;
realT maxabs= -REALmax;
realT sumabs= 0;
realT maxwidth= 0;
int k;
for (k= 0; k < dimension; k++) {
if (qh SCALElast && k == dimension-1)
abscoord= maxwidth;
else if (qh DELAUNAY && k == dimension-1) /* will qh_setdelaunay() */
abscoord= 2 * maxabs * maxabs; /* may be low by qh hull_dim/2 */
else {
maxcoord= -REALmax;
mincoord= REALmax;
FORALLpoint_(points, numpoints) {
maximize_(maxcoord, point[k]);
minimize_(mincoord, point[k]);
}
maximize_(maxwidth, maxcoord-mincoord);
abscoord= fmax_(maxcoord, -mincoord);
}
sumabs += abscoord;
maximize_(maxabs, abscoord);
} /* for k */
distround= qh_distround (qh hull_dim, maxabs, sumabs);
joggle= distround * qh_JOGGLEdefault;
maximize_(joggle, REALepsilon * qh_JOGGLEdefault);
trace2((qh ferr, "qh_detjoggle: joggle=%2.2g maxwidth=%2.2g\n", joggle, maxwidth));
return joggle;
} /* detjoggle */
/*-<a href="qh-c.htm#geom"
>-------------------------------</a><a name="detroundoff">-</a>
qh_detroundoff()
determine maximum roundoff errors from
REALepsilon, REALmax, REALmin, qh.hull_dim, qh.MAXabs_coord,
qh.MAXsumcoord, qh.MAXwidth, qh.MINdenom_1
accounts for qh.SETroundoff, qh.RANDOMdist, qh MERGEexact
qh.premerge_cos, qh.postmerge_cos, qh.premerge_centrum,
qh.postmerge_centrum, qh.MINoutside,
qh_RATIOnearinside, qh_COPLANARratio, qh_WIDEcoplanar
returns:
sets qh.DISTround, etc. (see below)
appends precision constants to qh.qhull_options
see:
qh_maxmin() for qh.NEARzero
design:
determine qh.DISTround for distance computations
determine minimum denominators for qh_divzero
determine qh.ANGLEround for angle computations
adjust qh.premerge_cos,... for roundoff error
determine qh.ONEmerge for maximum error due to a single merge
determine qh.NEARinside, qh.MAXcoplanar, qh.MINvisible,
qh.MINoutside, qh.WIDEfacet
initialize qh.max_vertex and qh.minvertex
*/
void qh_detroundoff (void) {
qh_option ("_max-width", NULL, &qh MAXwidth);
if (!qh SETroundoff) {
qh DISTround= qh_distround (qh hull_dim, qh MAXabs_coord, qh MAXsumcoord);
if (qh RANDOMdist)
qh DISTround += qh RANDOMfactor * qh MAXabs_coord;
qh_option ("Error-roundoff", NULL, &qh DISTround);
}
qh MINdenom_1= fmax_(1.0/REALmax, REALmin);
qh MINdenom= qh MINdenom_1 * qh MAXabs_coord;
qh MINdenom_1_2= sqrt (qh MINdenom_1 * qh hull_dim) ; /* if will be normalized */
qh MINdenom_2= qh MINdenom_1_2 * qh MAXabs_coord;
/* for inner product */
qh ANGLEround= 1.01 * qh hull_dim * REALepsilon;
if (qh RANDOMdist)
qh ANGLEround += qh RANDOMfactor;
if (qh premerge_cos < REALmax/2) {
qh premerge_cos -= qh ANGLEround;
if (qh RANDOMdist)
qh_option ("Angle-premerge-with-random", NULL, &qh premerge_cos);
}
if (qh postmerge_cos < REALmax/2) {
qh postmerge_cos -= qh ANGLEround;
if (qh RANDOMdist)
qh_option ("Angle-postmerge-with-random", NULL, &qh postmerge_cos);
}
qh premerge_centrum += 2 * qh DISTround; /*2 for centrum and distplane()*/
qh postmerge_centrum += 2 * qh DISTround;
if (qh RANDOMdist && (qh MERGEexact || qh PREmerge))
qh_option ("Centrum-premerge-with-random", NULL, &qh premerge_centrum);
if (qh RANDOMdist && qh POSTmerge)
qh_option ("Centrum-postmerge-with-random", NULL, &qh postmerge_centrum);
{ /* compute ONEmerge, max vertex offset for merging simplicial facets */
realT maxangle= 1.0, maxrho;
minimize_(maxangle, qh premerge_cos);
minimize_(maxangle, qh postmerge_cos);
/* max diameter * sin theta + DISTround for vertex to its hyperplane */
qh ONEmerge= sqrt (qh hull_dim) * qh MAXwidth *
sqrt (1.0 - maxangle * maxangle) + qh DISTround;
maxrho= qh hull_dim * qh premerge_centrum + qh DISTround;
maximize_(qh ONEmerge, maxrho);
maxrho= qh hull_dim * qh postmerge_centrum + qh DISTround;
maximize_(qh ONEmerge, maxrho);
if (qh MERGING)
qh_option ("_one-merge", NULL, &qh ONEmerge);
}
qh NEARinside= qh ONEmerge * qh_RATIOnearinside; /* only used if qh KEEPnearinside */
if (qh JOGGLEmax < REALmax/2 && (qh KEEPcoplanar || qh KEEPinside)) {
realT maxdist;
qh KEEPnearinside= True;
maxdist= sqrt (qh hull_dim) * qh JOGGLEmax + qh DISTround;
maxdist= 2*maxdist; /* vertex and coplanar point can joggle in opposite directions */
maximize_(qh NEARinside, maxdist); /* must agree with qh_nearcoplanar() */
}
if (qh KEEPnearinside)
qh_option ("_near-inside", NULL, &qh NEARinside);
if (qh JOGGLEmax < qh DISTround) {
ivp_message( "qhull error: the joggle for 'QJn', %.2g, is below roundoff for distance computations, %.2g\n",
qh JOGGLEmax, qh DISTround);
qh_errexit (qh_ERRinput, NULL, NULL);
}
if (qh MINvisible > REALmax/2) {
if (!qh MERGING)
qh MINvisible= qh DISTround;
else if (qh hull_dim <= 3)
qh MINvisible= qh premerge_centrum;
else
qh MINvisible= qh_COPLANARratio * qh premerge_centrum;
if (qh APPROXhull && qh MINvisible > qh MINoutside)
qh MINvisible= qh MINoutside;
qh_option ("Visible-distance", NULL, &qh MINvisible);
}
if (qh MAXcoplanar > REALmax/2) {
qh MAXcoplanar= qh MINvisible;
qh_option ("U-coplanar-distance", NULL, &qh MAXcoplanar);
}
if (!qh APPROXhull) { /* user may specify qh MINoutside */
qh MINoutside= 2 * qh MINvisible;
if (qh premerge_cos < REALmax/2)
maximize_(qh MINoutside, (1- qh premerge_cos) * qh MAXabs_coord);
qh_option ("Width-outside", NULL, &qh MINoutside);
}
qh WIDEfacet= qh MINoutside;
maximize_(qh WIDEfacet, qh_WIDEcoplanar * qh MAXcoplanar);
maximize_(qh WIDEfacet, qh_WIDEcoplanar * qh MINvisible);
qh_option ("_wide-facet", NULL, &qh WIDEfacet);
if (qh MINvisible > qh MINoutside + 3 * REALepsilon
&& !qh BESToutside && !qh FORCEoutput)
ivp_message( "qhull input warning: minimum visibility V%.2g is greater than \nminimum outside W%.2g. Flipped facets are likely.\n",
qh MINvisible, qh MINoutside);
qh max_vertex= qh DISTround;
qh min_vertex= -qh DISTround;
/* numeric constants reported in printsummary */
} /* detroundoff */
/*-<a href="qh-c.htm#geom"
>-------------------------------</a><a name="detsimplex">-</a>
qh_detsimplex( apex, points, dim, nearzero )
compute determinant of a simplex with point apex and base points
returns:
signed determinant and nearzero from qh_determinant
notes:
uses qh.gm_matrix/qh.gm_row (assumes they're big enough)
design:
construct qm_matrix by subtracting apex from points
compute determinate
*/
realT qh_detsimplex(pointT *apex, setT *points, int dim, boolT *nearzero) {
pointT *coorda, *coordp, *gmcoord, *point, **pointp;
coordT **rows;
int k, i=0;
realT det;
zinc_(Zdetsimplex);
gmcoord= qh gm_matrix;
rows= qh gm_row;
FOREACHpoint_(points) {
if (i == dim)
break;
rows[i++]= gmcoord;
coordp= point;
coorda= apex;
for (k= dim; k--; )
*(gmcoord++)= *coordp++ - *coorda++;
}
if (i < dim) {
ivp_message( "qhull internal error (qh_detsimplex): #points %d < dimension %d\n",
i, dim);
qh_errexit (qh_ERRqhull, NULL, NULL);
}
det= qh_determinant (rows, dim, nearzero);
trace2((qh ferr, "qh_detsimplex: det=%2.2g for point p%d, dim %d, nearzero? %d\n",
det, qh_pointid(apex), dim, *nearzero));
return det;
} /* detsimplex */
/*-<a href="qh-c.htm#geom"
>-------------------------------</a><a name="distnorm">-</a>
qh_distnorm( dim, point, normal, offset )
return distance from point to hyperplane at normal/offset
returns:
dist
notes:
dist > 0 if point is outside of hyperplane
see:
qh_distplane in geom.c
*/
realT qh_distnorm (int dim, pointT *point, pointT *normal, realT *offsetp) {
coordT *normalp= normal, *coordp= point;
realT dist;
int k;
dist= *offsetp;
for (k= dim; k--; )
dist += *(coordp++) * *(normalp++);
return dist;
} /* distnorm */
/*-<a href="qh-c.htm#geom"
>-------------------------------</a><a name="distround">-</a>
qh_distround ( dimension, maxabs, maxsumabs )
compute maximum round-off error for a distance computation
to a normalized hyperplane
maxabs is the maximum absolute value of a coordinate
maxsumabs is the maximum possible sum of absolute coordinate values
returns:
max dist round for REALepsilon
notes:
calculate roundoff error according to
Lemma 3.2-1 of Golub and van Loan "Matrix Computation"
use sqrt(dim) since one vector is normalized
or use maxsumabs since one vector is < 1
*/
realT qh_distround (int dimension, realT maxabs, realT maxsumabs) {
realT maxdistsum, maxround;
maxdistsum= sqrt (dimension) * maxabs;
minimize_( maxdistsum, maxsumabs);
maxround= REALepsilon * (dimension * maxdistsum * 1.01 + maxabs);
/* adds maxabs for offset */
trace4((qh ferr, "qh_distround: %2.2g maxabs %2.2g maxsumabs %2.2g maxdistsum %2.2g\n",
maxround, maxabs, maxsumabs, maxdistsum));
return maxround;
} /* distround */
/*-<a href="qh-c.htm#geom"
>-------------------------------</a><a name="divzero">-</a>
qh_divzero( numer, denom, mindenom1, zerodiv )
divide by a number that's nearly zero
mindenom1= minimum denominator for dividing into 1.0
returns:
quotient
sets zerodiv and returns 0.0 if it would overflow
design:
if numer is nearly zero and abs(numer) < abs(denom)
return numer/denom
else if numer is nearly zero
return 0 and zerodiv
else if denom/numer non-zero
return numer/denom
else
return 0 and zerodiv
*/
realT qh_divzero (realT numer, realT denom, realT mindenom1, boolT *zerodiv) {
realT temp, numerx, denomx;
if (numer < mindenom1 && numer > -mindenom1) {
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -