---------------------------------
geom2.c
infrequently used geometric routines of qhull
see qh-geom.htm and geom.h
copyright (c) 1993-2003 The Geometry Center
frequently used code goes into geom.c
*/
#include "qhull_a.h"
/*================== functions in alphabetic order ============*/
/*---------------------------------
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 */
/*---------------------------------
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 */
/*---------------------------------
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) {
fprintf (qh ferr, "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 */
/*---------------------------------
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 */
/*---------------------------------
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= 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; /* adjust qh.NEARinside for joggle */
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) {
fprintf (qh ferr, "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)
fprintf (qh ferr, "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 */
/*---------------------------------
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) {
fprintf (qh ferr, "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 */
/*---------------------------------
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 */
/*---------------------------------
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 */
/*---------------------------------
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) {
numerx= fabs_(numer);
denomx= fabs_(denom);
if (numerx < denomx) {
*zerodiv= False;
return numer/denom;
}else {
*zerodiv= True;
return 0.0;
}
}
temp= denom/numer;
if (temp > mindenom1 || temp < -mindenom1) {
*zerodiv= False;
return numer/denom;
}else {
*zerodiv= True;
return 0.0;
}
} /* divzero */
/*---------------------------------
qh_facetarea( facet )
return area for a facet
notes:
if non-simplicial,
uses centrum to triangulate facet and sums the projected areas.
if (qh DELAUNAY),
computes projected area instead for last coordinate
assumes facet->normal exists
projecting tricoplanar facets to the hyperplane does not appear to make a difference
design:
if simplicial
compute area
else
for each ridge
compute area from centrum to ridge
negate area if upper Delaunay facet
*/
realT qh_facetarea (facetT *facet) {
vertexT *apex;
pointT *centrum;
realT area= 0.0;
ridgeT *ridge, **ridgep;
if (facet->simplicial) {
apex= SETfirstt_(facet->vertices, vertexT);
area= qh_facetarea_simplex (qh hull_dim, apex->point, facet->vertices,
apex, facet->toporient, facet->normal, &facet->offset);
}else {
if (qh CENTERtype == qh_AScentrum)
centrum= facet->center;
else
centrum= qh_getcentrum (facet);
FOREACHridge_(facet->ridges)
area += qh_facetarea_simplex (qh hull_dim, centrum, ridge->vertices,
NULL, (ridge->top == facet), facet->normal, &facet->offset);
if (qh CENTERtype != qh_AScentrum)
qh_memfree (centrum, qh normal_size);
}
if (facet->upperdelaunay && qh DELAUNAY)
area= -area; /* the normal should be [0,...,1] */
trace4((qh ferr, "qh_facetarea: f%d area %2.2g\n", facet->id, area));
return area;
} /* facetarea */
/*---------------------------------
qh_facetarea_simplex( dim, apex, vertices, notvertex, toporient, normal, offset )
return area for a simplex defined by
an apex, a base of vertices, an orientation, and a unit normal
if simplicial or tricoplanar facet,
notvertex is defined and it is skipped in vertices
returns:
computes area of simplex projected to plane [normal,offset]
returns 0 if vertex too far below plane (qh WIDEfacet)
vertex can't be apex of tricoplanar facet
notes:
if (qh DELAUNAY),
computes projected area instead for last coordinate
uses qh gm_matrix/gm_row and qh hull_dim
helper function for qh_facetarea
design:
if Notvertex
translate simplex to apex
else
project simplex to normal/offset
translate simplex to apex
if Delaunay
set last row/column to 0 with -1 on diagonal
else
set last row to Normal
compute determinate
scale and flip sign for area
*/
realT qh_facetarea_simplex (int dim, coordT *apex, setT *vertices,
vertexT *notvertex, boolT toporient, coordT *normal, realT *offset) {
pointT *coorda, *coordp, *gmcoord;
coordT **rows, *normalp;
int k, i=0;
realT area, dist;
vertexT *vertex, **vertexp;
boolT nearzero;
gmcoord= qh gm_matrix;
rows= qh gm_row;
FOREACHvertex_(vertices) {
if (vertex == notvertex)
continue;
rows[i++]= gmcoord;
coorda= apex;
coordp= vertex->point;
normalp= normal;
if (notvertex) {
for (k= dim; k--; )
*(gmcoord++)= *coordp++ - *coorda++;
}else {
dist= *offset;
for (k= dim; k--; )
dist += *coordp++ * *normalp++;
if (dist < -qh WIDEfacet) {
zinc_(Znoarea);
return 0.0;
}
coordp= vertex->point;
normalp= normal;
for (k= dim; k--; )
*(gmcoord++)= (*coordp++ - dist * *normalp++) - *coorda++;
}
}
if (i != dim-1) {
fprintf (qh ferr, "qhull internal error (qh_facetarea_simplex): #points %d != dim %d -1\n",
i, dim);
qh_errexit (qh_ERRqhull, NULL, NULL);
}
rows[i]= gmcoord;
if (qh DELAUNAY) {
for (i= 0; i < dim-1; i++)
rows[i][dim-1]= 0.0;
for (k= dim; k--; )
*(gmcoord++)= 0.0;
rows[dim-1][dim-1]= -1.0;
}else {
normalp= normal;
for (k= dim; k--; )
*(gmcoord++)= *normalp++;
}
zinc_(Zdetsimplex);
area= qh_determinant (rows, dim, &nearzero);
if (toporient)
area= -area;
area *= qh AREAfactor;
trace4((qh ferr, "qh_facetarea_simplex: area=%2.2g for point p%d, toporient %d, nearzero? %d\n",
area, qh_pointid(apex), toporient, nearzero));
return area;
} /* facetarea_simplex */
/*---------------------------------
qh_facetcenter( vertices )
return Voronoi center (Voronoi vertex) for a facet's vertices
returns:
return temporary point equal to the center
see:
qh_voronoi_center()
*/
pointT *qh_facetcenter (setT *vertices) {
setT *points= qh_settemp (qh_setsize (vertices));
vertexT *vertex, **vertexp;
pointT *center;
FOREACHvertex_(vertices)
qh_setappend (&points, vertex->point);
center= qh_voronoi_center (qh hull_dim-1, points);
qh_settempfree (&points);
return center;
} /* facetcenter */
/*---------------------------------
qh_findgooddist( point, facetA, dist, facetlist )
find best good facet visible for point from facetA
assumes facetA is visible from point
returns:
best facet, i.e., good facet that is furthest from point
distance to best facet
NULL if none
moves good, visible facets (and some other visible facets)
to end of qh facet_list
notes:
uses qh visit_id
design:
initialize bestfacet if facetA is good
move facetA to end of facetlist
for each facet on facetlist
for each unvisited neighbor of facet
move visible neighbors to end of facetlist
update best good neighbor
if no good neighbors, update best facet
*/
facetT *qh_findgooddist (pointT *point, facetT *facetA, realT *distp,
facetT **facetlist) {
realT bestdist= -REALmax, dist;
facetT *neighbor, **neighborp, *bestfacet=NULL, *facet;
boolT goodseen= False;
if (facetA->good) {
zinc_(Zcheckpart); /* calls from check_bestdist occur after print stats */
qh_distplane (point, facetA, &bestdist);
bestfacet= facetA;
goodseen= True;
}
qh_removefacet (facetA);
qh_appendfacet (facetA);
*facetlist= facetA;
facetA->visitid= ++qh visit_id;
FORALLfacet_(*facetlist) {
FOREACHneighbor_(facet) {
if (neighbor->visitid == qh visit_id)
continue;
neighbor->visitid= qh visit_id;
if (goodseen && !neighbor->good)
continue;
zinc_(Zcheckpart);
qh_distplane (point, neighbor, &dist);
if (dist > 0) {
qh_removefacet (neighbor);
qh_appendfacet (neighbor);
if (neighbor->good) {
goodseen= True;
if (dist > bestdist) {
bestdist= dist;
bestfacet= neighbor;
}
}
}
}
}
if (bestfacet) {
*distp= bestdist;
trace2((qh ferr, "qh_findgooddist: p%d is %2.2g above good facet f%d\n",
qh_pointid(point), bestdist, bestfacet->id));
return bestfacet;
}
trace4((qh ferr, "qh_findgooddist: no good facet for p%d above f%d\n",
qh_pointid(point), facetA->id));
return NULL;
} /* findgooddist */
/*---------------------------------
qh_getarea( facetlist )
set area of all facets in facetlist
collect statistics
returns:
sets qh totarea/totvol to total area and volume of convex hull
for Delaunay triangulation, computes projected area of the lower or upper hull
ignores upper hull if qh ATinfinity
notes:
could compute outer volume by expanding facet area by rays from interior
the following attempt at perpendicular projection underestimated badly:
qh.totoutvol += (-dist + facet->maxoutside + qh DISTround)
* area/ qh hull_dim;
design:
for each facet on facetlist
compute facet->area
update qh.totarea and qh.totvol
*/
void qh_getarea (facetT *facetlist) {
realT area;
realT dist;
facetT *facet;
if (qh REPORTfreq)
fprintf (qh ferr, "computing area of each facet and volume of the convex hull\n");
else
trace1((qh ferr, "qh_getarea: computing volume and area for each facet\n"));
qh totarea= qh totvol= 0.0;
FORALLfacet_(facetlist) {
if (!facet->normal)
continue;
if (facet->upperdelaunay && qh ATinfinity)
continue;
facet->f.area= area= qh_facetarea (facet);
facet->isarea= True;
if (qh DELAUNAY) {
if (facet->upperdelaunay == qh UPPERdelaunay)
qh totarea += area;
}else {
qh totarea += area;
qh_distplane (qh interior_point, facet, &dist);
qh totvol += -dist * area/ qh hull_dim;
}
if (qh PRINTstatistics) {
wadd_(Wareatot, area);
wmax_(Wareamax, area);
wmin_(Wareamin, area);
}
}
} /* getarea */
/*---------------------------------
qh_gram_schmidt( dim, row )
implements Gram-Schmidt orthogonalization by rows
returns:
false if zero norm
overwrites rows[dim][dim]
notes:
see Golub & van Loan Algorithm 6.2-2
overflow due to small divisors not handled
design:
for each row
compute norm for row
if non-zero, normalize row
for each remaining rowA
compute inner product of row and rowA
reduce rowA by row * inner product
*/
boolT qh_gram_schmidt(int dim, realT **row) {
realT *rowi, *rowj, norm;
int i, j, k;
for(i=0; i < dim; i++) {
rowi= row[i];
for (norm= 0.0, k= dim; k--; rowi++)
norm += *rowi * *rowi;
norm= sqrt(norm);
wmin_(Wmindenom, norm);
if (norm == 0.0) /* either 0 or overflow due to sqrt */
return False;
for(k= dim; k--; )
*(--rowi) /= norm;
for(j= i+1; j < dim; j++) {
rowj= row[j];
for(norm= 0.0, k=dim; k--; )
norm += *rowi++ * *rowj++;
for(k=dim; k--; )
*(--rowj) -= *(--rowi) * norm;
}
}
return True;
} /* gram_schmidt */
/*---------------------------------
qh_inthresholds( normal, angle )
return True if normal within qh.lower_/upper_threshold
returns:
estimate of angle by summing of threshold diffs
angle may be NULL
smaller "angle" is better
notes:
invalid if qh.SPLITthresholds
see:
qh.lower_threshold in qh_initbuild()
qh_initthresholds()
design:
for each dimension
test threshold
*/
boolT qh_inthresholds (coordT *normal, realT *angle) {
boolT within= True;
int k;
realT threshold;
if (angle)
*angle= 0.0;
for(k= 0; k < qh hull_dim; k++) {
threshold= qh lower_threshold[k];
if (threshold > -REALmax/2) {
if (normal[k] < threshold)
within= False;
if (angle) {
threshold -= normal[k];
*angle += fabs_(threshold);
}
}
if (qh upper_threshold[k] < REALmax/2) {
threshold= qh upper_threshold[k];
if (normal[k] > threshold)
within= False;
if (angle) {
threshold -= normal[k];
*angle += fabs_(threshold);
}
}
}
return within;
} /* inthresholds */
/*---------------------------------
qh_joggleinput()
randomly joggle input to Qhull by qh.JOGGLEmax
initial input is qh.first_point/qh.num_points of qh.hull_dim
repeated calls use qh.input_points/qh.num_points
returns:
joggles points at qh.first_point/qh.num_points
copies data to qh.input_points/qh.input_malloc if first time
determines qh.JOGGLEmax if it was zero
if qh.DELAUNAY
computes the Delaunay projection of the joggled points
notes:
if qh.DELAUNAY, unnecessarily joggles the last coordinate
the initial 'QJn' may be set larger than qh_JOGGLEmaxincrease
design:
if qh.DELAUNAY
set qh.SCALElast for reduced precision errors
if first call
initialize qh.input_points to the original input points
if qh.JOGGLEmax == 0
determine default qh.JOGGLEmax
else
increase qh.JOGGLEmax according to qh.build_cnt
joggle the input by adding a random number in [-qh.JOGGLEmax,qh.JOGGLEmax]
if qh.DELAUNAY
sets the Delaunay projection
*/
void qh_joggleinput (void) {
int size, i, seed;
coordT *coordp, *inputp;
realT randr, randa, randb;
if (!qh input_points) { /* first call */
qh input_points= qh first_point;
qh input_malloc= qh POINTSmalloc;
size= qh num_points * qh hull_dim * sizeof(coordT);
if (!(qh first_point=(coordT*)malloc(size))) {
fprintf(qh ferr, "qhull error: insufficient memory to joggle %d points\n",
qh num_points);
qh_errexit(qh_ERRmem, NULL, NULL);
}
qh POINTSmalloc= True;
if (qh JOGGLEmax == 0.0) {
qh JOGGLEmax= qh_detjoggle (qh input_points, qh num_points, qh hull_dim);
qh_option ("QJoggle", NULL, &qh JOGGLEmax);
}
}else { /* repeated call */
if (!qh RERUN && qh build_cnt > qh_JOGGLEretry) {
if (((qh build_cnt-qh_JOGGLEretry-1) % qh_JOGGLEagain) == 0) {
realT maxjoggle= qh MAXwidth * qh_JOGGLEmaxincrease;
if (qh JOGGLEmax < maxjoggle) {
qh JOGGLEmax *= qh_JOGGLEincrease;
minimize_(qh JOGGLEmax, maxjoggle);
}
}
}
qh_option ("QJoggle", NULL, &qh JOGGLEmax);
}
if (qh build_cnt > 1 && qh JOGGLEmax > fmax_(qh MAXwidth/4, 0.1)) {
fprintf (qh ferr, "qhull error: the current joggle for 'QJn', %.2g, is too large for the width\nof the input. If possible, recompile Qhull with higher-precision reals.\n",
qh JOGGLEmax);
qh_errexit (qh_ERRqhull, NULL, NULL);
}
/* for some reason, using qh ROTATErandom and qh_RANDOMseed does not repeat the run. Use 'TRn' instead */
seed= qh_RANDOMint;
qh_option ("_joggle-seed", &seed, NULL);
trace0((qh ferr, "qh_joggleinput: joggle input by %2.2g with seed %d\n",
qh JOGGLEmax, seed));
inputp= qh input_points;
coordp= qh first_point;
randa= 2.0 * qh JOGGLEmax/qh_RANDOMmax;
randb= -qh JOGGLEmax;
size= qh num_points * qh hull_dim;
for (i= size; i--; ) {
randr= qh_RANDOMint;
*(coordp++)= *(inputp++) + (randr * randa + randb);
}
if (qh DELAUNAY) {
qh last_low= qh last_high= qh last_newhigh= REALmax;
qh_setdelaunay (qh hull_dim, qh num_points, qh first_point);
}
} /* joggleinput */
/*---------------------------------
qh_maxabsval( normal, dim )
return pointer to maximum absolute value of a dim vector
returns NULL if dim=0
*/
realT *qh_maxabsval (realT *normal, int dim) {
realT maxval= -REALmax;
realT *maxp= NULL, *colp, absval;
int k;
for (k= dim, colp= normal; k--; colp++) {
absval= fabs_(*colp);
if (absval > maxval) {
maxval= absval;
maxp= colp;
}
}
return maxp;
} /* maxabsval */
/*---------------------------------
qh_maxmin( points, numpoints, dimension )
return max/min points for each dimension
determine max and min coordinates
returns:
returns a temporary set of max and min points
may include duplicate points. Does not include qh.GOO}t
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) {
numerx= fabs_(numer);
denomx= fabs_(denom);
if (numerx < denomx) {
*zerodiv= False;
return numer/denom;
}else {
*zerodiv= True;
return 0.0;
}
}
temp= denom/numer;
if (temp > mindenom1 || temp < -mindenom1) {
*zerodiv= False;
return numer/denom;
}else {
*zerodiv= True;
return 0.0;
}
} /* divzero */
/*---------------------------------
qh_facetarea( facet )
return area for a facet
notes:
if non-simplicial,
uses centrum to triangulate facet and sums the projected areas.
if (qh DELAUNAY),
computes projected area instead for last coordinate
assumes facet->normal exists
projecting tricoplanar facets to the hyperplane does not appear to make a difference
design:
if simplicial
compute area
else
for each ridge
compute area from centrum to ridge
negate area if upper Delaunay facet
*/
realT qh_facetarea (facetT *facet) {
vertexT *apex;
pointT *centrum;
realT area= 0.0;
ridgeT *ridge, **ridgep;
if (facet->simplicial) {
apex= SETfirstt_(facet->vertices, vertexT);
area= qh_facetarea_simplex (qh hull_dim, apex->point, facet->vertices,
apex, facet->toporient, facet->normal, &facet->offset);
}else {
if (qh CENTERtype == qh_AScentrum)
centrum= facet->center;
else
centrum= qh_getcentrum (facet);
FOREACHridge_(facet->ridges)
area += qh_facetarea_simplex (qh hull_dim, centrum, ridge->vertices,
NULL, (ridge->top == facet), facet->normal, &facet->offset);
if (qh CENTERtype != qh_AScentrum)
qh_memfree (centrum, qh normal_size);
}
if (facet->upperdelaunay && qh DELAUNAY)
area= -area; /* the normal should be [0,...,1] */
trace4((qh ferr, "qh_facetarea: f%d area %2.2g\n", facet->id, area));
return area;
} /* facetarea */
/*---------------------------------
qh_facetarea_simplex( dim, apex, vertices, notvertex, toporient, normal, offset )
return area for a simplex defined by
an apex, a base of vertices, an orientation, and a unit normal
if simplicial or tricoplanar facet,
notvertex is defined and it is skipped in vertices
returns:
computes area of simplex projected to plane [normal,offset]
returns 0 if vertex too far below plane (qh WIDEfacet)
vertex can't be apex of tricoplanar facet
notes:
if (qh DELAUNAY),
computes projected area instead for last coordinate
uses qh gm_matrix/gm_row and qh hull_dim
helper function for qh_facetarea
design:
if Notvertex
translate simplex to apex
else
project simplex to normal/offset
translate simplex to apex
if Delaunay
set last row/column to 0 with -1 on diagonal
else
set last row to Normal
compute determinate
scale and flip sign for area
*/
realT qh_facetarea_simplex (int dim, coordT *apex, setT *vertices,
vertexT *notvertex, boolT toporient, coordT *normal, realT *offset) {
pointT *coorda, *coordp, *gmcoord;
coordT **rows, *normalp;
int k, i=0;
realT area, dist;
vertexT *vertex, **vertexp;
boolT nearzero;
gmcoord= qh gm_matrix;
rows= qh gm_row;
FOREACHvertex_(vertices) {
if (vertex == notvertex)
continue;
rows[i++]= gmcoord;
coorda= apex;
coordp= vertex->point;
normalp= normal;
if (notvertex) {
for (k= dim; k--; )
*(gmcoord++)= *coordp++ - *coorda++;
}else {
dist= *offset;
for (k= dim; k--; )
dist += *coordp++ * *normalp++;
if (dist < -qh WIDEfacet) {
zinc_(Znoarea);
return 0.0;
}
coordp= vertex->point;
normalp= normal;
for (k= dim; k--; )
*(gmcoord++)= (*coordp++ - dist * *normalp++) - *coorda++;
}
}
if (i != dim-1) {
fprintf (qh ferr, "qhull internal error (qh_facetarea_simplex): #points %d != dim %d -1\n",
i, dim);
qh_errexit (qh_ERRqhull, NULL, NULL);
}
rows[i]= gmcoord;
if (qh DELAUNAY) {
for (i= 0; i < dim-1; i++)
rows[i][dim-1]= 0.0;
for (k= dim; k--; )
*(gmcoord++)= 0.0;
rows[dim-1][dim-1]= -1.0;
}else {
normalp= normal;
for (k= dim; k--; )
*(gmcoord++)= *normalp++;
}
zinc_(Zdetsimplex);
area= qh_determinant (rows, dim, &nearzero);
if (toporient)
area= -area;
area *= qh AREAfactor;
trace4((qh ferr, "qh_facetarea_simplex: area=%2.2g for point p%d, toporient %d, nearzero? %d\n",
area, qh_pointid(apex), toporient, nearzero));
return area;
} /* facetarea_simplex */
/*---------------------------------
qh_facetcenter( vertices )
return Voronoi center (Voronoi vertex) for a facet's vertices
returns:
return temporary point equal to the center
see:
qh_voronoi_center()
*/
pointT *qh_facetcenter (setT *vertices) {
setT *points= qh_settemp (qh_setsize (vertices));
vertexT *vertex, **vertexp;
pointT *center;
FOREACHvertex_(vertices)
qh_setappend (&points, vertex->point);
center= qh_voronoi_center (qh hull_dim-1, points);
qh_settempfree (&points);
return center;
} /* facetcenter */
/*---------------------------------
qh_findgooddist( point, facetA, dist, facetlist )
find best good facet visible for point from facetA
assumes facetA is visible from point
returns:
best facet, i.e., good facet that is furthest from point
distance to best facet
NULL if none
moves good, visible facets (and some other visible facets)
to end of qh facet_list
notes:
uses qh visit_id
design:
initialize bestfacet if facetA is good
move facetA to end of facetlist
for each facet on facetlist
for each unvisited neighbor of facet
move visible neighbors to end of facetlist
update best good neighbor
if no good neighbors, update best facet
*/
facetT *qh_findgooddist (pointT *point, facetT *facetA, realT *distp,
facetT **facetlist) {
realT bestdist= -REALmax, dist;
facetT *neighbor, **neighborp, *bestfacet=NULL, *facet;
boolT goodseen= False;
if (facetA->good) {
zinc_(Zcheckpart); /* calls from check_bestdist occur after print stats */
qh_distplane (point, facetA, &bestdist);
bestfacet= facetA;
goodseen= True;
}
qh_removefacet (facetA);
qh_appendfacet (facetA);
*facetlist= facetA;
facetA->visitid= ++qh visit_id;
FORALLfacet_(*facetlist) {
FOREACHneighbor_(facet) {
if (neighbor->visitid == qh visit_id)
continue;
neighbor->visitid= qh visit_id;
if (goodseen && !neighbor->good)
continue;
zinc_(Zcheckpart);
qh_distplane (point, neighbor, &dist);
if (dist > 0) {
qh_removefacet (neighbor);
qh_appendfacet (neighbor);
if (neighbor->good) {
goodseen= True;
if (dist > bestdist) {
bestdist= dist;
bestfacet= neighbor;
}
}
}
}
}
if (bestfacet) {
*distp= bestdist;
trace2((qh ferr, "qh_findgooddist: p%d is %2.2g above good facet f%d\n",
qh_pointid(point), bestdist, bestfacet->id));
return bestfacet;
}
trace4((qh ferr, "qh_findgooddist: no good facet for p%d above f%d\n",
qh_pointid(point), facetA->id));
return NULL;
} /* findgooddist */
/*---------------------------------
qh_getarea( facetlist )
set area of all facets in facetlist
collect statistics
returns:
sets qh totarea/totvol to total area and volume of convex hull
for Delaunay triangulation, computes projected area of the lower or upper hull
ignores upper hull if qh ATinfinity
notes:
could compute outer volume by expanding facet area by rays from interior
the following attempt at perpendicular projection underestimated badly:
qh.totoutvol += (-dist + facet->maxoutside + qh DISTround)
* area/ qh hull_dim;
design:
for each facet on facetlist
compute facet->area
update qh.totarea and qh.totvol
*/
void qh_getarea (facetT *facetlist) {
realT area;
realT dist;
facetT *facet;
if (qh REPORTfreq)
fprintf (qh ferr, "computing area of each facet and volume of the convex hull\n");
else
trace1((qh ferr, "qh_getarea: computing volume and area for each facet\n"));
qh totarea= qh totvol= 0.0;
FORALLfacet_(facetlist) {
if (!facet->normal)
continue;
if (facet->upperdelaunay && qh ATinfinity)
continue;
facet->f.area= area= qh_facetarea (facet);
facet->isarea= True;
if (qh DELAUNAY) {
if (facet->upperdelaunay == qh UPPERdelaunay)
qh totarea += area;
}else {
qh totarea += area;
qh_distplane (qh interior_point, facet, &dist);
qh totvol += -dist * area/ qh hull_dim;
}
if (qh PRINTstatistics) {
wadd_(Wareatot, area);
wmax_(Wareamax, area);
wmin_(Wareamin, area);
}
}
} /* getarea */
/*---------------------------------
qh_gram_schmidt( dim, row )
implements Gram-Schmidt orthogonalization by rows
returns:
false if zero norm
overwrites rows[dim][dim]
notes:
see Golub & van Loan Algorithm 6.2-2
overflow due to small divisors not handled
design:
for each row
compute norm for row
if non-zero, normalize row
for each remaining rowA
compute inner product of row and rowA
reduce rowA by row * inner product
*/
boolT qh_gram_schmidt(int dim, realT **row) {
realT *rowi, *rowj, norm;
int i, j, k;
for(i=0; i < dim; i++) {
rowi= row[i];
for (norm= 0.0, k= dim; k--; rowi++)
norm += *rowi * *rowi;
norm= sqrt(norm);
wmin_(Wmindenom, norm);
if (norm == 0.0) /* either 0 or overflow due to sqrt */
return False;
for(k= dim; k--; )
*(--rowi) /= norm;
for(j= i+1; j < dim; j++) {
rowj= row[j];
for(norm= 0.0, k=dim; k--; )
norm += *rowi++ * *rowj++;
for(k=dim; k--; )
*(--rowj) -= *(--rowi) * norm;
}
}
return True;
} /* gram_schmidt */
/*---------------------------------
qh_inthresholds( normal, angle )
return True if normal within qh.lower_/upper_threshold
returns:
estimate of angle by summing of threshold diffs
angle may be NULL
smaller "angle" is better
notes:
invalid if qh.SPLITthresholds
see:
qh.lower_threshold in qh_initbuild()
qh_initthresholds()
design:
for each dimension
test threshold
*/
boolT qh_inthresholds (coordT *normal, realT *angle) {
boolT within= True;
int k;
realT threshold;
if (angle)
*angle= 0.0;
for(k= 0; k < qh hull_dim; k++) {
threshold= qh lower_threshold[k];
if (threshold > -REALmax/2) {
if (normal[k] < threshold)
within= False;
if (angle) {
threshold -= normal[k];
*angle += fabs_(threshold);
}
}
if (qh upper_threshold[k] < REALmax/2) {
threshold= qh upper_threshold[k];
if (normal[k] > threshold)
within= False;
if (angle) {
threshold -= normal[k];
*angle += fabs_(threshold);
}
}
}
return within;
} /* inthresholds */
/*---------------------------------
qh_joggleinput()
randomly joggle input to Qhull by qh.JOGGLEmax
initial input is qh.first_point/qh.num_points of qh.hull_dim
repeated calls use qh.input_points/qh.num_points
returns:
joggles points at qh.first_point/qh.num_points
copies data to qh.input_points/qh.input_malloc if first time
determines qh.JOGGLEmax if it was zero
if qh.DELAUNAY
computes the Delaunay projection of the joggled points
notes:
if qh.DELAUNAY, unnecessarily joggles the last coordinate
the initial 'QJn' may be set larger than qh_JOGGLEmaxincrease
design:
if qh.DELAUNAY
set qh.SCALElast for reduced precision errors
if first call
initialize qh.input_points to the original input points
if qh.JOGGLEmax == 0
determine default qh.JOGGLEmax
else
increase qh.JOGGLEmax according to qh.build_cnt
joggle the input by adding a random number in [-qh.JOGGLEmax,qh.JOGGLEmax]
if qh.DELAUNAY
sets the Delaunay projection
*/
void qh_joggleinput (void) {
int size, i, seed;
coordT *coordp, *inputp;
realT randr, randa, randb;
if (!qh input_points) { /* first call */
qh input_points= qh first_point;
qh input_malloc= qh POINTSmalloc;
size= qh num_points * qh hull_dim * sizeof(coordT);
if (!(qh first_point=(coordT*)malloc(size))) {
fprintf(qh ferr, "qhull error: insufficient memory to joggle %d points\n",
qh num_points);
qh_errexit(qh_ERRmem, NULL, NULL);
}
qh POINTSmalloc= True;
if (qh JOGGLEmax == 0.0) {
qh JOGGLEmax= qh_detjoggle (qh input_points, qh num_points, qh hull_dim);
qh_option ("QJoggle", NULL, &qh JOGGLEmax);
}
}else { /* repeated call */
if (!qh RERUN && qh build_cnt > qh_JOGGLEretry) {
if (((qh build_cnt-qh_JOGGLEretry-1) % qh_JOGGLEagain) == 0) {
realT maxjoggle= qh MAXwidth * qh_JOGGLEmaxincrease;
if (qh JOGGLEmax < maxjoggle) {
qh JOGGLEmax *= qh_JOGGLEincrease;
minimize_(qh JOGGLEmax, maxjoggle);
}
}
}
qh_option ("QJoggle", NULL, &qh JOGGLEmax);
}
if (qh build_cnt > 1 && qh JOGGLEmax > fmax_(qh MAXwidth/4, 0.1)) {
fprintf (qh ferr, "qhull error: the current joggle for 'QJn', %.2g, is too large for the width\nof the input. If possible, recompile Qhull with higher-precision reals.\n",
qh JOGGLEmax);
qh_errexit (qh_ERRqhull, NULL, NULL);
}
/* for some reason, using qh ROTATErandom and qh_RANDOMseed does not repeat the run. Use 'TRn' instead */
seed= qh_RANDOMint;
qh_option ("_joggle-seed", &seed, NULL);
trace0((qh ferr, "qh_joggleinput: joggle input by %2.2g with seed %d\n",
qh JOGGLEmax, seed));
inputp= qh input_points;
coordp= qh first_point;
randa= 2.0 * qh JOGGLEmax/qh_RANDOMmax;
randb= -qh JOGGLEmax;
size= qh num_points * qh hull_dim;
for (i= size; i--; ) {
randr= qh_RANDOMint;
*(coordp++)= *(inputp++) + (randr * randa + randb);
}
if (qh DELAUNAY) {
qh last_low= qh last_high= qh last_newhigh= REALmax;
qh_setdelaunay (qh hull_dim, qh num_points, qh first_point);
}
} /* joggleinput */
/*---------------------------------
qh_maxabsval( normal, dim )
return pointer to maximum absolute value of a dim vector
returns NULL if dim=0
*/
realT *qh_maxabsval (realT *normal, int dim) {
realT maxval= -REALmax;
realT *maxp= NULL, *colp, absval;
int k;
for (k= dim, colp= normal; k--; colp++) {
absval= fabs_(*colp);
if (absval > maxval) {
maxval= absval;
maxp= colp;
}
}
return maxp;
} /* maxabsval */
/*---------------------------------
qh_maxmin( points, numpoints, dimension )
return max/min points for each dimension
determine max and min coordinates
returns:
returns a temporary set of max and min points
may include duplicate points. Does not include qh.GOO}t
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) {
numerx= fabs_(numer);
denomx= fabs_(denom);
if (numerx < denomx) {
*zerodiv= False;
return numer/denom;
}else {
*zerodiv= True;
return 0.0;
}
}
temp= denom/numer;
if (temp > mindenom1 || temp < -mindenom1) {
*zerodiv= False;
return numer/denom;
}else {
*zerodiv= True;
return 0.0;
}
} /* divzero */
/*---------------------------------
qh_facetarea( facet )
return area for a facet
notes:
if non-simplicial,
uses centrum to triangulate facet and sums the projected areas.
if (qh DELAUNAY),
computes projected area instead for last coordinate
assumes facet->normal exists
projecting tricoplanar facets to the hyperplane does not appear to make a difference
design:
if simplicial
compute area
else
for each ridge
compute area from centrum to ridge
negate area if upper Delaunay facet
*/
realT qh_facetarea (facetT *facet) {
vertexT *apex;
pointT *centrum;
realT area= 0.0;
ridgeT *ridge, **ridgep;
if (facet->simplicial) {
apex= SETfirstt_(facet->vertices, vertexT);
area= qh_facetarea_simplex (qh hull_dim, apex->point, facet->vertices,
apex, facet->toporient, facet->normal, &facet->offset);
}else {
if (qh CENTERtype == qh_AScentrum)
centrum= facet->center;
else
centrum= qh_getcentrum (facet);
FOREACHridge_(facet->ridges)
area += qh_facetarea_simplex (qh hull_dim, centrum, ridge->vertices,
NULL, (ridge->top == facet), facet->normal, &facet->offset);
if (qh CENTERtype != qh_AScentrum)
qh_memfree (centrum, qh normal_size);
}
if (facet->upperdelaunay && qh DELAUNAY)
area= -area; /* the normal should be [0,...,1] */
trace4((qh ferr, "qh_facetarea: f%d area %2.2g\n", facet->id, area));
return area;
} /* facetarea */
/*---------------------------------
qh_facetarea_simplex( dim, apex, vertices, notvertex, toporient, normal, offset )
return area for a simplex defined by
an apex, a base of vertices, an orientation, and a unit normal
if simplicial or tricoplanar facet,
notvertex is defined and it is skipped in vertices
returns:
computes area of simplex projected to plane [normal,offset]
returns 0 if vertex too far below plane (qh WIDEfacet)
vertex can't be apex of tricoplanar facet
notes:
if (qh DELAUNAY),
computes projected area instead for last coordinate
uses qh gm_matrix/gm_row and qh hull_dim
helper function for qh_facetarea
design:
if Notvertex
translate simplex to apex
else
project simplex to normal/offset
translate simplex to apex
if Delaunay
set last row/column to 0 with -1 on diagonal
else
set last row to Normal
compute determinate
scale and flip sign for area
*/
realT qh_facetarea_simplex (int dim, coordT *apex, setT *vertices,
vertexT *notvertex, boolT toporient, coordT *normal, realT *offset) {
pointT *coorda, *coordp, *gmcoord;
coordT **rows, *normalp;
int k, i=0;
realT area, dist;
vertexT *vertex, **vertexp;
boolT nearzero;
gmcoord= qh gm_matrix;
rows= qh gm_row;
FOREACHvertex_(vertices) {
if (vertex == notvertex)
continue;
rows[i++]= gmcoord;
coorda= apex;
coordp= vertex->point;
normalp= normal;
if (notvertex) {
for (k= dim; k--; )
*(gmcoord++)= *coordp++ - *coorda++;
}else {
dist= *offset;
for (k= dim; k--; )
dist += *coordp++ * *normalp++;
if (dist < -qh WIDEfacet) {
zinc_(Znoarea);
return 0.0;
}
coordp= vertex->point;
normalp= normal;
for (k= dim; k--; )
*(gmcoord++)= (*coordp++ - dist * *normalp++) - *coorda++;
}
}
if (i != dim-1) {
fprintf (qh ferr, "qhull internal error (qh_facetarea_simplex): #points %d != dim %d -1\n",
i, dim);
qh_errexit (qh_ERRqhull, NULL, NULL);
}
rows[i]= gmcoord;
if (qh DELAUNAY) {
for (i= 0; i < dim-1; i++)
rows[i][dim-1]= 0.0;
for (k= dim; k--; )
*(gmcoord++)= 0.0;
rows[dim-1][dim-1]= -1.0;
}else {
normalp= normal;
for (k= dim; k--; )
*(gmcoord++)= *normalp++;
}
zinc_(Zdetsimplex);
area= qh_determinant (rows, dim, &nearzero);
if (toporient)
area= -area;
area *= qh AREAfactor;
trace4((qh ferr, "qh_facetarea_simplex: area=%2.2g for point p%d, toporient %d, nearzero? %d\n",
area, qh_pointid(apex), toporient, nearzero));
return area;
} /* facetarea_simplex */
/*---------------------------------
qh_facetcenter( vertices )
return Voronoi center (Voronoi vertex) for a facet's vertices
returns:
return temporary point equal to the center
see:
qh_voronoi_center()
*/
pointT *qh_facetcenter (setT *vertices) {
setT *points= qh_settemp (qh_setsize (vertices));
vertexT *vertex, **vertexp;
pointT *center;
FOREACHvertex_(vertices)
qh_setappend (&points, vertex->point);
center= qh_voronoi_center (qh hull_dim-1, points);
qh_settempfree (&points);
return center;
} /* facetcenter */
/*---------------------------------
qh_findgooddist( point, facetA, dist, facetlist )
find best good facet visible for point from facetA
assumes facetA is visible from point
returns:
best facet, i.e., good facet that is furthest from point
distance to best facet
NULL if none
moves good, visible facets (and some other visible facets)
to end of qh facet_list
notes:
uses qh visit_id
design:
initialize bestfacet if facetA is good
move facetA to end of facetlist
for each facet on facetlist
for each unvisited neighbor of facet
move visible neighbors to end of facetlist
update best good neighbor
if no good neighbors, update best facet
*/
facetT *qh_findgooddist (pointT *point, facetT *facetA, realT *distp,
facetT **facetlist) {
realT bestdist= -REALmax, dist;
facetT *neighbor, **neighborp, *bestfacet=NULL, *facet;
boolT goodseen= False;
if (facetA->good) {
zinc_(Zcheckpart); /* calls from check_bestdist occur after print stats */
qh_distplane (point, facetA, &bestdist);
bestfacet= facetA;
goodseen= True;
}
qh_removefacet (facetA);
qh_appendfacet (facetA);
*facetlist= facetA;
facetA->visitid= ++qh visit_id;
FORALLfacet_(*facetlist) {
FOREACHneighbor_(facet) {
if (neighbor->visitid == qh visit_id)
continue;
neighbor->visitid= qh visit_id;
if (goodseen && !neighbor->good)
continue;
zinc_(Zcheckpart);
qh_distplane (point, neighbor, &dist);
if (dist > 0) {
qh_removefacet (neighbor);
qh_appendfacet (neighbor);
if (neighbor->good) {
goodseen= True;
if (dist > bestdist) {
bestdist= dist;
bestfacet= neighbor;
}
}
}
}
}
if (bestfacet) {
*distp= bestdist;
trace2((qh ferr, "qh_findgooddist: p%d is %2.2g above good facet f%d\n",
qh_pointid(point), bestdist, bestfacet->id));
return bestfacet;
}
trace4((qh ferr, "qh_findgooddist: no good facet for p%d above f%d\n",
qh_pointid(point), facetA->id));
return NULL;
} /* findgooddist */
/*---------------------------------
qh_getarea( facetlist )
set area of all facets in facetlist
collect statistics
returns:
sets qh totarea/totvol to total area and volume of convex hull
for Delaunay triangulation, computes projected area of the lower or upper hull
ignores upper hull if qh ATinfinity
notes:
could compute outer volume by expanding facet area by rays from interior
the following attempt at perpendicular projection underestimated badly:
qh.totoutvol += (-dist + facet->maxoutside + qh DISTround)
* area/ qh hull_dim;
design:
for each facet on facetlist
compute facet->area
update qh.totarea and qh.totvol
*/
void qh_getarea (facetT *facetlist) {
realT area;
realT dist;
facetT *facet;
if (qh REPORTfreq)
fprintf (qh ferr, "computing area of each facet and volume of the convex hull\n");
else
trace1((qh ferr, "qh_getarea: computing volume and area for each facet\n"));
qh totarea= qh totvol= 0.0;
FORALLfacet_(facetlist) {
if (!facet->normal)
continue;
if (facet->upperdelaunay && qh ATinfinity)
continue;
facet->f.area= area= qh_facetarea (facet);
facet->isarea= True;
if (qh DELAUNAY) {
if (facet->upperdelaunay == qh UPPERdelaunay)
qh totarea += area;
}else {
qh totarea += area;
qh_distplane (qh interior_point, facet, &dist);
qh totvol += -dist * area/ qh hull_dim;
}
if (qh PRINTstatistics) {
wadd_(Wareatot, area);
wmax_(Wareamax, area);
wmin_(Wareamin, area);
}
}
} /* getarea */
/*---------------------------------
qh_gram_schmidt( dim, row )
implements Gram-Schmidt orthogonalization by rows
returns:
false if zero norm
overwrites rows[dim][dim]
notes:
see Golub & van Loan Algorithm 6.2-2
overflow due to small divisors not handled
design:
for each row
compute norm for row
if non-zero, normalize row
for each remaining rowA
compute inner product of row and rowA
reduce rowA by row * inner product
*/
boolT qh_gram_schmidt(int dim, realT **row) {
realT *rowi, *rowj, norm;
int i, j, k;
for(i=0; i < dim; i++) {
rowi= row[i];
for (norm= 0.0, k= dim; k--; rowi++)
norm += *rowi * *rowi;
norm= sqrt(norm);
wmin_(Wmindenom, norm);
if (norm == 0.0) /* either 0 or overflow due to sqrt */
return False;
for(k= dim; k--; )
*(--rowi) /= norm;
for(j= i+1; j < dim; j++) {
rowj= row[j];
for(norm= 0.0, k=dim; k--; )
norm += *rowi++ * *rowj++;
for(k=dim; k--; )
*(--rowj) -= *(--rowi) * norm;
}
}
return True;
} /* gram_schmidt */
/*---------------------------------
qh_inthresholds( normal, angle )
return True if normal within qh.lower_/upper_threshold
returns:
estimate of angle by summing of threshold diffs
angle may be NULL
smaller "angle" is better
notes:
invalid if qh.SPLITthresholds
see:
qh.lower_threshold in qh_initbuild()
qh_initthresholds()
design:
for each dimension
test threshold
*/
boolT qh_inthresholds (coordT *normal, realT *angle) {
boolT within= True;
int k;
realT threshold;
if (angle)
*angle= 0.0;
for(k= 0; k < qh hull_dim; k++) {
threshold= qh lower_threshold[k];
if (threshold > -REALmax/2) {
if (normal[k] < threshold)
within= False;
if (angle) {
threshold -= normal[k];
*angle += fabs_(threshold);
}
}
if (qh upper_threshold[k] < REALmax/2) {
threshold= qh upper_threshold[k];
if (normal[k] > threshold)
within= False;
if (angle) {
threshold -= normal[k];
*angle += fabs_(threshold);
}
}
}
return within;
} /* inthresholds */
/*---------------------------------
qh_joggleinput()
randomly joggle input to Qhull by qh.JOGGLEmax
initial input is qh.first_point/qh.num_points of qh.hull_dim
repeated calls use qh.input_points/qh.num_points
returns:
joggles points at qh.first_point/qh.num_points
copies data to qh.input_points/qh.input_malloc if first time
determines qh.JOGGLEmax if it was zero
if qh.DELAUNAY
computes the Delaunay projection of the joggled points
notes:
if qh.DELAUNAY, unnecessarily joggles the last coordinate
the initial 'QJn' may be set larger than qh_JOGGLEmaxincrease
design:
if qh.DELAUNAY
set qh.SCALElast for reduced precision errors
if first call
initialize qh.input_points to the original input points
if qh.JOGGLEmax == 0
determine default qh.JOGGLEmax
else
increase qh.JOGGLEmax according to qh.build_cnt
joggle the input by adding a random number in [-qh.JOGGLEmax,qh.JOGGLEmax]
if qh.DELAUNAY
sets the Delaunay projection
*/
void qh_joggleinput (void) {
int size, i, seed;
coordT *coordp, *inputp;
realT randr, randa, randb;
if (!qh input_points) { /* first call */
qh input_points= qh first_point;
qh input_malloc= qh POINTSmalloc;
size= qh num_points * qh hull_dim * sizeof(coordT);
if (!(qh first_point=(coordT*)malloc(size))) {
fprintf(qh ferr, "qhull error: insufficient memory to joggle %d points\n",
qh num_points);
qh_errexit(qh_ERRmem, NULL, NULL);
}
qh POINTSmalloc= True;
if (qh JOGGLEmax == 0.0) {
qh JOGGLEmax= qh_detjoggle (qh input_points, qh num_points, qh hull_dim);
qh_option ("QJoggle", NULL, &qh JOGGLEmax);
}
}else { /* repeated call */
if (!qh RERUN && qh build_cnt > qh_JOGGLEretry) {
if (((qh build_cnt-qh_JOGGLEretry-1) % qh_JOGGLEagain) == 0) {
realT maxjoggle= qh MAXwidth * qh_JOGGLEmaxincrease;
if (qh JOGGLEmax < maxjoggle) {
qh JOGGLEmax *= qh_JOGGLEincrease;
minimize_(qh JOGGLEmax, maxjoggle);
}
}
}
qh_option ("QJoggle", NULL, &qh JOGGLEmax);
}
if (qh build_cnt > 1 && qh JOGGLEmax > fmax_(qh MAXwidth/4, 0.1)) {
fprintf (qh ferr, "qhull error: the current joggle for 'QJn', %.2g, is too large for the width\nof the input. If possible, recompile Qhull with higher-precision reals.\n",
qh JOGGLEmax);
qh_errexit (qh_ERRqhull, NULL, NULL);
}
/* for some reason, using qh ROTATErandom and qh_RANDOMseed does not repeat the run. Use 'TRn' instead */
seed= qh_RANDOMint;
qh_option ("_joggle-seed", &seed, NULL);
trace0((qh ferr, "qh_joggleinput: joggle input by %2.2g with seed %d\n",
qh JOGGLEmax, seed));
inputp= qh input_points;
coordp= qh first_point;
randa= 2.0 * qh JOGGLEmax/qh_RANDOMmax;
randb= -qh JOGGLEmax;
size= qh num_points * qh hull_dim;
for (i= size; i--; ) {
randr= qh_RANDOMint;
*(coordp++)= *(inputp++) + (randr * randa + randb);
}
if (qh DELAUNAY) {
qh last_low= qh last_high= qh last_newhigh= REALmax;
qh_setdelaunay (qh hull_dim, qh num_points, qh first_point);
}
} /* joggleinput */
/*---------------------------------
qh_maxabsval( normal, dim )
return pointer to maximum absolute value of a dim vector
returns NULL if dim=0
*/
realT *qh_maxabsval (realT *normal, int dim) {
realT maxval= -REALmax;
realT *maxp= NULL, *colp, absval;
int k;
for (k= dim, colp= normal; k--; colp++) {
absval= fabs_(*colp);
if (absval > maxval) {
maxval= absval;
maxp= colp;
}
}
return maxp;
} /* maxabsval */
/*---------------------------------
qh_maxmin( points, numpoints, dimension )
return max/min points for each dimension
determine max and min coordinates
returns:
returns a temporary set of max and min points
may include duplicate points. Does not include qh.GOO}t
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) {
numerx= fabs_(numer);
denomx= fabs_(denom);
if (numerx < denomx) {
*zerodiv= False;
return numer/denom;
}else {
*zerodiv= True;
return 0.0;
}
}
temp= denom/numer;
if (temp > mindenom1 || temp < -mindenom1) {
*zerodiv= False;
return numer/denom;
}else {
*zerodiv= True;
return 0.0;
}
} /* divzero */
/*---------------------------------
qh_facetarea( facet )
return area for a facet
notes:
if non-simplicial,
uses centrum to triangulate facet and sums the projected areas.
if (qh DELAUNAY),
computes projected area instead for last coordinate
assumes facet->normal exists
projecting tricoplanar facets to the hyperplane does not appear to make a difference
design:
if simplicial
compute area
else
for each ridge
compute area from centrum to ridge
negate area if upper Delaunay facet
*/
realT qh_facetarea (facetT *facet) {
vertexT *apex;
pointT *centrum;
realT area= 0.0;
ridgeT *ridge, **ridgep;
if (facet->simplicial) {
apex= SETfirstt_(facet->vertices, vertexT);
area= qh_facetarea_simplex (qh hull_dim, apex->point, facet->vertices,
apex, facet->toporient, facet->normal, &facet->offset);
}else {
if (qh CENTERtype == qh_AScentrum)
centrum= facet->center;
else
centrum= qh_getcentrum (facet);
FOREACHridge_(facet->ridges)
area += qh_facetarea_simplex (qh hull_dim, centrum, ridge->vertices,
NULL, (ridge->top == facet), facet->normal, &facet->offset);
if (qh CENTERtype != qh_AScentrum)
qh_memfree (centrum, qh normal_size);
}
if (facet->upperdelaunay && qh DELAUNAY)
area= -area; /* the normal should be [0,...,1] */
trace4((qh ferr, "qh_facetarea: f%d area %2.2g\n", facet->id, area));
return area;
} /* facetarea */
/*---------------------------------
qh_facetarea_simplex( dim, apex, vertices, notvertex, toporient, normal, offset )
return area for a simplex defined by
an apex, a base of vertices, an orientation, and a unit normal
if simplicial or tricoplanar facet,
notvertex is defined and it is skipped in vertices
returns:
computes area of simplex projected to plane [normal,offset]
returns 0 if vertex too far below plane (qh WIDEfacet)
vertex can't be apex of tricoplanar facet
notes:
if (qh DELAUNAY),
computes projected area instead for last coordinate
uses qh gm_matrix/gm_row and qh hull_dim
helper function for qh_facetarea
design:
if Notvertex
translate simplex to apex
else
project simplex to normal/offset
translate simplex to apex
if Delaunay
set last row/column to 0 with -1 on diagonal
else
set last row to Normal
compute determinate
scale and flip sign for area
*/
realT qh_facetarea_simplex (int dim, coordT *apex, setT *vertices,
vertexT *notvertex, boolT toporient, coordT *normal, realT *offset) {
pointT *coorda, *coordp, *gmcoord;
coordT **rows, *normalp;
int k, i=0;
realT area, dist;
vertexT *vertex, **vertexp;
boolT nearzero;
gmcoord= qh gm_matrix;
rows= qh gm_row;
FOREACHvertex_(vertices) {
if (vertex == notvertex)
continue;
rows[i++]= gmcoord;
coorda= apex;
coordp= vertex->point;
normalp= normal;
if (notvertex) {
for (k= dim; k--; )
*(gmcoord++)= *coordp++ - *coorda++;
}else {
dist= *offset;
for (k= dim; k--; )
dist