/*
--------------------------------- geom2.c infrequently used geometric routines of qhull see qh-geom.htm and geom.h Copyright (c) 1993-2020 The Geometry Center. $Id: //main/2019/qhull/src/libqhull/geom2.c#14 $$Change: 3037 $ $DateTime: 2020/09/03 17:28:32 $$Author: bbarber $ frequently used code goes into geom.c */ #include "qhull_a.h" /*================== functions in alphabetic order ============*/ /*--------------------------------- qh_copypoints( points, numpoints, dimension ) return qh_malloc'd copy of points notes: qh_free the returned points to avoid a memory leak */ coordT *qh_copypoints(coordT *points, int numpoints, int dimension) { int size; coordT *newpoints; size= numpoints * dimension * (int)sizeof(coordT); if (!(newpoints= (coordT *)qh_malloc((size_t)size))) { qh_fprintf(qh ferr, 6004, "qhull error: insufficient memory to copy %d points\n", numpoints); qh_errexit(qh_ERRmem, NULL, NULL); } memcpy((char *)newpoints, (char *)points, (size_t)size); /* newpoints!=0 by QH6004 */ 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] (!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) { qh_fprintf(qh ferr, 6005, "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) < 10*qh NEARzero[1]) /* QH11031 FIX: 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) < 10*qh NEARzero[2]) /* QH11031 FIX: what should this be? det 5.5e-12 was flat for qh_maxsimplex of qdelaunay 0,0 27,27 -36,36 -9,63 */ *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; if (qh SETroundoff) distround= qh DISTround; /* 'En' */ else{ 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, 2001, "qh_detjoggle: joggle=%2.2g maxwidth=%2.2g\n", joggle, maxwidth)); return joggle; } /* detjoggle */ /*--------------------------------- qh_detmaxoutside(); determine qh.MAXoutside target for qh_RATIO... tests of distance updates option '_max-outside' notes: called from qh_addpoint and qh_detroundoff accounts for qh.ONEmerge, qh.DISTround, qh.MINoutside ('Wn'), qh.max_outside see qh_maxout for qh.max_outside with qh.DISTround */ void qh_detmaxoutside(void) { realT maxoutside; maxoutside= fmax_(qh max_outside, qh ONEmerge + qh DISTround); maximize_(maxoutside, qh MINoutside); qh MAXoutside= maxoutside; trace3((qh ferr, 3056, "qh_detmaxoutside: MAXoutside %2.2g from qh.max_outside %2.2g, ONEmerge %2.2g, MINoutside %2.2g, DISTround %2.2g\n", qh MAXoutside, qh max_outside, qh ONEmerge, qh MINoutside, qh DISTround)); } /* detmaxoutside */ /*--------------------------------- 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); 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; trace4((qh ferr, 4096, "qh_detroundoff: increase qh.ANGLEround by option 'R%2.2g'\n", 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((realT)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((realT)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) { qh_fprintf(qh ferr, 6006, "qhull option 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-max-coplanar", 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) qh_fprintf(qh ferr, 7001, "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 */ qh_detmaxoutside(); } /* 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: called by qh_maxsimplex and qh_initialvertices 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) { qh_fprintf(qh ferr, 6007, "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, 2002, "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 if qh.RANDOMdist ('Qr'), adjusts qh_distround returns: max dist round for qh.REALepsilon and qh.RANDOMdist notes: calculate roundoff error according to Golub & van Loan, 1983, Lemma 3.2-1, "Rounding Errors" 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, delta; maxdistsum= sqrt((realT)dimension) * maxabs; minimize_( maxdistsum, maxsumabs); maxround= REALepsilon * (dimension * maxdistsum * 1.01 + maxabs); /* adds maxabs for offset */ if (qh RANDOMdist) { delta= qh RANDOMfactor * maxabs; maxround += delta; trace4((qh ferr, 4092, "qh_distround: increase roundoff by random delta %2.2g for option 'R%2.2g'\n", delta, qh RANDOMfactor)); } trace4((qh ferr, 4008, "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, (boolT)(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, 4009, "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) { qh_fprintf(qh ferr, 6008, "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_(Zdetfacetarea); area= qh_determinant(rows, dim, &nearzero); if (toporient) area= -area; area *= qh AREAfactor; trace4((qh ferr, 4010, "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) { zzinc_(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; zzinc_(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, 2003, "qh_findgooddist: p%d is %2.2g above good facet f%d\n", qh_pointid(point), bestdist, bestfacet->id)); return bestfacet; } trace4((qh ferr, 4011, "qh_findgooddist: no good facet for p%d above f%d\n", qh_pointid(point), facetA->id)); return NULL; } /* findgooddist */ /*--------------------------------- qh_furthestnewvertex( unvisited, facet, &maxdist ) return furthest unvisited, new vertex to a facet return: NULL if no vertex is above facet maxdist to facet updates v.visitid notes: Ignores vertices in facetB Does not change qh.vertex_visit. Use in conjunction with qh_furthestvertex */ vertexT *qh_furthestnewvertex(unsigned int unvisited, facetT *facet, realT *maxdistp /* qh.newvertex_list */) { vertexT *maxvertex= NULL, *vertex; coordT dist, maxdist= 0.0; FORALLvertex_(qh newvertex_list) { if (vertex->newfacet && vertex->visitid <= unvisited) { vertex->visitid= qh vertex_visit; qh_distplane(vertex->point, facet, &dist); if (dist > maxdist) { maxdist= dist; maxvertex= vertex; } } } trace4((qh ferr, 4085, "qh_furthestnewvertex: v%d dist %2.2g is furthest new vertex for f%d\n", getid_(maxvertex), maxdist, facet->id)); *maxdistp= maxdist; return maxvertex; } /* furthestnewvertex */ /*--------------------------------- qh_furthestvertex( facetA, facetB, &maxdist, &mindist ) return furthest vertex in facetA from facetB, or NULL if none return: maxdist and mindist to facetB or 0.0 if none updates qh.vertex_visit notes: Ignores vertices in facetB */ vertexT *qh_furthestvertex(facetT *facetA, facetT *facetB, realT *maxdistp, realT *mindistp) { vertexT *maxvertex= NULL, *vertex, **vertexp; coordT dist, maxdist= -REALmax, mindist= REALmax; qh vertex_visit++; FOREACHvertex_(facetB->vertices) vertex->visitid= qh vertex_visit; FOREACHvertex_(facetA->vertices) { if (vertex->visitid != qh vertex_visit) { vertex->visitid= qh vertex_visit; zzinc_(Zvertextests); qh_distplane(vertex->point, facetB, &dist); if (!maxvertex) { maxdist= dist; mindist= dist; maxvertex= vertex; }else if (dist > maxdist) { maxdist= dist; maxvertex= vertex; }else if (dist < mindist) mindist= dist; } } if (!maxvertex) { trace3((qh ferr, 3067, "qh_furthestvertex: all vertices of f%d are in f%d. Returning 0.0 for max and mindist\n", facetA->id, facetB->id)); maxdist= mindist= 0.0; }else { trace4((qh ferr, 4084, "qh_furthestvertex: v%d dist %2.2g is furthest (mindist %2.2g) of f%d above f%d\n", maxvertex->id, maxdist, mindist, facetA->id, facetB->id)); } *maxdistp= maxdist; *mindistp= mindist; return maxvertex; } /* furthestvertex */ /*--------------------------------- qh_getarea( facetlist ) set area of all facets in facetlist collect statistics nop if hasAreaVolume 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 hasAreaVolume) return; if (qh REPORTfreq) qh_fprintf(qh ferr, 8020, "computing area of each facet and volume of the convex hull\n"); else trace1((qh ferr, 1001, "qh_getarea: computing area for each facet and its volume to qh.interior_point (dist*area/dim)\n")); qh totarea= qh totvol= 0.0; FORALLfacet_(facetlist) { if (!facet->normal) continue; if (facet->upperdelaunay && qh ATinfinity) continue; if (!facet->isarea) { facet->f.area= qh_facetarea(facet); facet->isarea= True; } area= facet->f.area; 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); } } qh hasAreaVolume= True; } /* 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, 1983, Algorithm 6.2-2, "Modified Gram-Schmidt" 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 i, seed, size; 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 * (int)sizeof(coordT); if (!(qh first_point= (coordT *)qh_malloc((size_t)size))) { qh_fprintf(qh ferr, 6009, "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)) { qh_fprintf(qh ferr, 6010, "qhull input error (qh_joggleinput): 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_ERRinput, 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, 6, "qh_joggleinput: joggle input by %4.4g 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.GOODpoint sets qh.NEARzero, qh.MAXabs_coord, qh.MAXsumcoord, qh.MAXwidth qh.MAXlastcoord, qh.MINlastcoord initializes qh.max_outside, qh.min_vertex, qh.WAScoplanar, qh.ZEROall_ok notes: loop duplicated in qh_detjoggle() design: initialize global precision variables checks definition of REAL... for each dimension for each point collect maximum and minimum point collect maximum of maximums and minimum of minimums determine qh.NEARzero for Gaussian Elimination */ setT *qh_maxmin(pointT *points, int numpoints, int dimension) { int k; realT maxcoord, temp; pointT *minimum, *maximum, *point, *pointtemp; setT *set; qh max_outside= 0.0; qh MAXabs_coord= 0.0; qh MAXwidth= -REALmax; qh MAXsumcoord= 0.0; qh min_vertex= 0.0; qh WAScoplanar= False; if (qh ZEROcentrum) qh ZEROall_ok= True; if (REALmin < REALepsilon && REALmin < REALmax && REALmin > -REALmax && REALmax > 0.0 && -REALmax < 0.0) ; /* all ok */ else { qh_fprintf(qh ferr, 6011, "qhull error: one or more floating point constants in user.h are inconsistent. REALmin %g, -REALmax %g, 0.0, REALepsilon %g, REALmax %g\n", REALmin, -REALmax, REALepsilon, REALmax); qh_errexit(qh_ERRinput, NULL, NULL); } set= qh_settemp(2*dimension); trace1((qh ferr, 8082, "qh_maxmin: dim min max width nearzero min-point max-point\n")); for (k=0; k < dimension; k++) { if (points == qh GOODpointp) minimum= maximum= points + dimension; else minimum= maximum= points; FORALLpoint_(points, numpoints) { if (point == qh GOODpointp) continue; if (maximum[k] < point[k]) maximum= point; else if (minimum[k] > point[k]) minimum= point; } if (k == dimension-1) { qh MINlastcoord= minimum[k]; qh MAXlastcoord= maximum[k]; } if (qh SCALElast && k == dimension-1) maxcoord= qh MAXabs_coord; else { maxcoord= fmax_(maximum[k], -minimum[k]); if (qh GOODpointp) { temp= fmax_(qh GOODpointp[k], -qh GOODpointp[k]); maximize_(maxcoord, temp); } temp= maximum[k] - minimum[k]; maximize_(qh MAXwidth, temp); } maximize_(qh MAXabs_coord, maxcoord); qh MAXsumcoord += maxcoord; qh_setappend(&set, minimum); qh_setappend(&set, maximum); /* calculation of qh NEARzero is based on Golub & van Loan, 1983, Eq. 4.4-13 for "Gaussian elimination with complete pivoting". Golub & van Loan say that n^3 can be ignored and 10 be used in place of rho */ qh NEARzero[k]= 80 * qh MAXsumcoord * REALepsilon; trace1((qh ferr, 8106, " %3d % 14.8e % 14.8e % 14.8e %4.4e p%-9d p%-d\n", k, minimum[k], maximum[k], maximum[k]-minimum[k], qh NEARzero[k], qh_pointid(minimum), qh_pointid(maximum))); if (qh SCALElast && k == dimension-1) trace1((qh ferr, 8107, " last coordinate scaled to (%4.4g, %4.4g), width %4.4e for option 'Qbb'\n", qh MAXabs_coord - qh MAXwidth, qh MAXabs_coord, qh MAXwidth)); } if (qh IStracing >= 1) qh_printpoints(qh ferr, "qh_maxmin: found the max and min points (by dim):", set); return(set); } /* maxmin */ /*--------------------------------- qh_maxouter( ) return maximum distance from facet to outer plane normally this is qh.max_outside+qh.DISTround does not include qh.JOGGLEmax see: qh_outerinner() notes: need to add another qh.DISTround if testing actual point with computation see qh_detmaxoutside for a qh_RATIO... target for joggle: qh_setfacetplane() updated qh.max_outer for Wnewvertexmax (max distance to vertex) need to use Wnewvertexmax since could have a coplanar point for a high facet that is replaced by a low facet need to add qh.JOGGLEmax if testing input points */ realT qh_maxouter(void) { realT dist; dist= fmax_(qh max_outside, qh DISTround); dist += qh DISTround; trace4((qh ferr, 4012, "qh_maxouter: max distance from facet to outer plane is %4.4g, qh.max_outside is %4.4g\n", dist, qh max_outside)); return dist; } /* maxouter */ /*--------------------------------- qh_maxsimplex( dim, maxpoints, points, numpoints, simplex ) determines maximum simplex for a set of points maxpoints is the subset of points with a min or max coordinate may start with points already in simplex skips qh.GOODpointp (assumes that it isn't in maxpoints) returns: simplex with dim+1 points notes: called by qh_initialvertices, qh_detvnorm, and qh_voronoi_center requires qh.MAXwidth to estimate determinate for each vertex assumes at least needed points in points maximizes determinate for x,y,z,w, etc. uses maxpoints as long as determinate is clearly non-zero design: initialize simplex with at least two points (find points with max or min x coordinate) create a simplex of dim+1 vertices as follows add point from maxpoints that maximizes the determinate of the point and the simplex vertices if last point and maxdet/prevdet < qh_RATIOmaxsimplex (3.0e-2) flag maybe_falsenarrow if no maxpoint or maxnearzero or maybe_falsenarrow search all points for maximum determinate early exit if maybe_falsenarrow and !maxnearzero and maxdet > prevdet */ void qh_maxsimplex(int dim, setT *maxpoints, pointT *points, int numpoints, setT **simplex) { pointT *point, **pointp, *pointtemp, *maxpoint, *minx=NULL, *maxx=NULL; boolT nearzero, maxnearzero= False, maybe_falsenarrow; int i, sizinit; realT maxdet= -1.0, prevdet= -1.0, det, mincoord= REALmax, maxcoord= -REALmax, mindet, ratio, targetdet; if (qh MAXwidth <= 0.0) { qh_fprintf(qh ferr, 6421, "qhull internal error (qh_maxsimplex): qh.MAXwidth required for qh_maxsimplex. Used to estimate determinate\n"); qh_errexit(qh_ERRqhull, NULL, NULL); } sizinit= qh_setsize(*simplex); if (sizinit >= 2) { maxdet= pow(qh MAXwidth, sizinit - 1); }else { if (qh_setsize(maxpoints) >= 2) { FOREACHpoint_(maxpoints) { if (maxcoord < point[0]) { maxcoord= point[0]; maxx= point; } if (mincoord > point[0]) { mincoord= point[0]; minx= point; } } }else { FORALLpoint_(points, numpoints) { if (point == qh GOODpointp) continue; if (maxcoord < point[0]) { maxcoord= point[0]; maxx= point; } if (mincoord > point[0]) { mincoord= point[0]; minx= point; } } } maxdet= maxcoord - mincoord; qh_setunique(simplex, minx); if (qh_setsize(*simplex) < 2) qh_setunique(simplex, maxx); sizinit= qh_setsize(*simplex); if (sizinit < 2) { qh_joggle_restart("input has same x coordinate"); if (zzval_(Zsetplane) > qh hull_dim+1) { qh_fprintf(qh ferr, 6012, "qhull precision error (qh_maxsimplex for voronoi_center): %d points with the same x coordinate %4.4g\n", qh_setsize(maxpoints)+numpoints, mincoord); qh_errexit(qh_ERRprec, NULL, NULL); }else { qh_fprintf(qh ferr, 6013, "qhull input error: input is less than %d-dimensional since all points have the same x coordinate %4.4g\n", qh hull_dim, mincoord); qh_errexit(qh_ERRinput, NULL, NULL); } } } for (i=sizinit; i < dim+1; i++) { prevdet= maxdet; maxpoint= NULL; maxdet= -1.0; FOREACHpoint_(maxpoints) { if (!qh_setin(*simplex, point) && point != maxpoint) { det= qh_detsimplex(point, *simplex, i, &nearzero); /* retests maxpoints if duplicate or multiple iterations */ if ((det= fabs_(det)) > maxdet) { maxdet= det; maxpoint= point; maxnearzero= nearzero; } } } maybe_falsenarrow= False; ratio= 1.0; targetdet= prevdet * qh MAXwidth; mindet= 10 * qh_RATIOmaxsimplex * targetdet; if (maxdet > 0.0) { ratio= maxdet / targetdet; if (ratio < qh_RATIOmaxsimplex) maybe_falsenarrow= True; } if (!maxpoint || maxnearzero || maybe_falsenarrow) { zinc_(Zsearchpoints); if (!maxpoint) { trace0((qh ferr, 7, "qh_maxsimplex: searching all points for %d-th initial vertex, better than mindet %4.4g, targetdet %4.4g\n", i+1, mindet, targetdet)); }else if (qh ALLpoints) { trace0((qh ferr, 30, "qh_maxsimplex: searching all points ('Qs') for %d-th initial vertex, better than p%d det %4.4g, targetdet %4.4g, ratio %4.4g\n", i+1, qh_pointid(maxpoint), maxdet, targetdet, ratio)); }else if (maybe_falsenarrow) { trace0((qh ferr, 17, "qh_maxsimplex: searching all points for %d-th initial vertex, better than p%d det %4.4g and mindet %4.4g, ratio %4.4g\n", i+1, qh_pointid(maxpoint), maxdet, mindet, ratio)); }else { trace0((qh ferr, 8, "qh_maxsimplex: searching all points for %d-th initial vertex, better than p%d det %2.2g and mindet %4.4g, targetdet %4.4g\n", i+1, qh_pointid(maxpoint), maxdet, mindet, targetdet)); } FORALLpoint_(points, numpoints) { if (point == qh GOODpointp) continue; if (!qh_setin(maxpoints, point) && !qh_setin(*simplex, point)) { det= qh_detsimplex(point, *simplex, i, &nearzero); if ((det= fabs_(det)) > maxdet) { maxdet= det; maxpoint= point; maxnearzero= nearzero; if (!maxnearzero && !qh ALLpoints && maxdet > mindet) break; } } } } /* !maxpoint */ if (!maxpoint) { qh_fprintf(qh ferr, 6014, "qhull internal error (qh_maxsimplex): not enough points available\n"); qh_errexit(qh_ERRqhull, NULL, NULL); } qh_setappend(simplex, maxpoint); trace1((qh ferr, 1002, "qh_maxsimplex: selected point p%d for %d`th initial vertex, det=%4.4g, targetdet=%4.4g, mindet=%4.4g\n", qh_pointid(maxpoint), i+1, maxdet, prevdet * qh MAXwidth, mindet)); } /* i */ } /* maxsimplex */ /*--------------------------------- qh_minabsval( normal, dim ) return minimum absolute value of a dim vector */ realT qh_minabsval(realT *normal, int dim) { realT minval= 0; realT maxval= 0; realT *colp; int k; for (k=dim, colp=normal; k--; colp++) { maximize_(maxval, *colp); minimize_(minval, *colp); } return fmax_(maxval, -minval); } /* minabsval */ /*--------------------------------- qh_mindif( vecA, vecB, dim ) return index of min abs. difference of two vectors */ int qh_mindiff(realT *vecA, realT *vecB, int dim) { realT mindiff= REALmax, diff; realT *vecAp= vecA, *vecBp= vecB; int k, mink= 0; for (k=0; k < dim; k++) { diff= *vecAp++ - *vecBp++; diff= fabs_(diff); if (diff < mindiff) { mindiff= diff; mink= k; } } return mink; } /* mindiff */ /*--------------------------------- qh_orientoutside( facet ) make facet outside oriented via qh.interior_point returns: True if facet reversed orientation. */ boolT qh_orientoutside(facetT *facet) { int k; realT dist; qh_distplane(qh interior_point, facet, &dist); if (dist > 0) { for (k=qh hull_dim; k--; ) facet->normal[k]= -facet->normal[k]; facet->offset= -facet->offset; return True; } return False; } /* orientoutside */ /*--------------------------------- qh_outerinner( facet, outerplane, innerplane ) if facet and qh.maxoutdone (i.e., qh_check_maxout) returns outer and inner plane for facet else returns maximum outer and inner plane accounts for qh.JOGGLEmax see: qh_maxouter(), qh_check_bestdist(), qh_check_points() notes: outerplaner or innerplane may be NULL facet is const Does not error (QhullFacet) includes qh.DISTround for actual points adds another qh.DISTround if testing with floating point arithmetic */ void qh_outerinner(facetT *facet, realT *outerplane, realT *innerplane) { realT dist, mindist; vertexT *vertex, **vertexp; if (outerplane) { if (!qh_MAXoutside || !facet || !qh maxoutdone) { *outerplane= qh_maxouter(); /* includes qh.DISTround */ }else { /* qh_MAXoutside ... */ #if qh_MAXoutside *outerplane= facet->maxoutside + qh DISTround; #endif } if (qh JOGGLEmax < REALmax/2) *outerplane += qh JOGGLEmax * sqrt((realT)qh hull_dim); } if (innerplane) { if (facet) { mindist= REALmax; FOREACHvertex_(facet->vertices) { zinc_(Zdistio); qh_distplane(vertex->point, facet, &dist); minimize_(mindist, dist); } *innerplane= mindist - qh DISTround; }else *innerplane= qh min_vertex - qh DISTround; if (qh JOGGLEmax < REALmax/2) *innerplane -= qh JOGGLEmax * sqrt((realT)qh hull_dim); } } /* outerinner */ /*--------------------------------- qh_pointdist( point1, point2, dim ) return distance between two points notes: returns distance squared if 'dim' is negative */ coordT qh_pointdist(pointT *point1, pointT *point2, int dim) { coordT dist, diff; int k; dist= 0.0; for (k= (dim > 0 ? dim : -dim); k--; ) { diff= *point1++ - *point2++; dist += diff * diff; } if (dim > 0) return(sqrt(dist)); return dist; } /* pointdist */ /*--------------------------------- qh_printmatrix( fp, string, rows, numrow, numcol ) print matrix to fp given by row vectors print string as header notes: print a vector by qh_printmatrix(fp, "", &vect, 1, len) */ void qh_printmatrix(FILE *fp, const char *string, realT **rows, int numrow, int numcol) { realT *rowp; realT r; /*bug fix*/ int i,k; qh_fprintf(fp, 9001, "%s\n", string); for (i=0; i < numrow; i++) { rowp= rows[i]; for (k=0; k < numcol; k++) { r= *rowp++; qh_fprintf(fp, 9002, "%6.3g ", r); } qh_fprintf(fp, 9003, "\n"); } } /* printmatrix */ /*--------------------------------- qh_printpoints( fp, string, points ) print pointids to fp for a set of points if string, prints string and 'p' point ids */ void qh_printpoints(FILE *fp, const char *string, setT *points) { pointT *point, **pointp; if (string) { qh_fprintf(fp, 9004, "%s", string); FOREACHpoint_(points) qh_fprintf(fp, 9005, " p%d", qh_pointid(point)); qh_fprintf(fp, 9006, "\n"); }else { FOREACHpoint_(points) qh_fprintf(fp, 9007, " %d", qh_pointid(point)); qh_fprintf(fp, 9008, "\n"); } } /* printpoints */ /*--------------------------------- qh_projectinput( ) project input points using qh.lower_bound/upper_bound and qh DELAUNAY if qh.lower_bound[k]=qh.upper_bound[k]= 0, removes dimension k if halfspace intersection removes dimension k from qh.feasible_point input points in qh first_point, num_points, input_dim returns: new point array in qh first_point of qh hull_dim coordinates sets qh POINTSmalloc if qh DELAUNAY projects points to paraboloid lowbound/highbound is also projected if qh ATinfinity adds point "at-infinity" if qh POINTSmalloc frees old point array notes: checks that qh.hull_dim agrees with qh.input_dim, PROJECTinput, and DELAUNAY design: sets project[k] to -1 (delete), 0 (keep), 1 (add for Delaunay) determines newdim and newnum for qh hull_dim and qh num_points projects points to newpoints projects qh.lower_bound to itself projects qh.upper_bound to itself if qh DELAUNAY if qh ATINFINITY projects points to paraboloid computes "infinity" point as vertex average and 10% above all points else uses qh_setdelaunay to project points to paraboloid */ void qh_projectinput(void) { int k,i; int newdim= qh input_dim, newnum= qh num_points; signed char *project; int projectsize= (qh input_dim + 1) * (int)sizeof(*project); pointT *newpoints, *coord, *infinity; realT paraboloid, maxboloid= 0; project= (signed char *)qh_memalloc(projectsize); memset((char *)project, 0, (size_t)projectsize); for (k=0; k < qh input_dim; k++) { /* skip Delaunay bound */ if (qh lower_bound[k] == 0.0 && qh upper_bound[k] == 0.0) { project[k]= -1; newdim--; } } if (qh DELAUNAY) { project[k]= 1; newdim++; if (qh ATinfinity) newnum++; } if (newdim != qh hull_dim) { qh_memfree(project, projectsize); qh_fprintf(qh ferr, 6015, "qhull internal error (qh_projectinput): dimension after projection %d != hull_dim %d\n", newdim, qh hull_dim); qh_errexit(qh_ERRqhull, NULL, NULL); } if (!(newpoints= qh temp_malloc= (coordT *)qh_malloc((size_t)(newnum * newdim) * sizeof(coordT)))) { qh_memfree(project, projectsize); qh_fprintf(qh ferr, 6016, "qhull error: insufficient memory to project %d points\n", qh num_points); qh_errexit(qh_ERRmem, NULL, NULL); } /* qh_projectpoints throws error if mismatched dimensions */ qh_projectpoints(project, qh input_dim+1, qh first_point, qh num_points, qh input_dim, newpoints, newdim); trace1((qh ferr, 1003, "qh_projectinput: updating lower and upper_bound\n")); qh_projectpoints(project, qh input_dim+1, qh lower_bound, 1, qh input_dim+1, qh lower_bound, newdim+1); qh_projectpoints(project, qh input_dim+1, qh upper_bound, 1, qh input_dim+1, qh upper_bound, newdim+1); if (qh HALFspace) { if (!qh feasible_point) { qh_memfree(project, projectsize); qh_fprintf(qh ferr, 6017, "qhull internal error (qh_projectinput): HALFspace defined without qh.feasible_point\n"); qh_errexit(qh_ERRqhull, NULL, NULL); } qh_projectpoints(project, qh input_dim, qh feasible_point, 1, qh input_dim, qh feasible_point, newdim); } qh_memfree(project, projectsize); if (qh POINTSmalloc) qh_free(qh first_point); qh first_point= newpoints; qh POINTSmalloc= True; qh temp_malloc= NULL; if (qh DELAUNAY && qh ATinfinity) { coord= qh first_point; infinity= qh first_point + qh hull_dim * qh num_points; for (k=qh hull_dim-1; k--; ) infinity[k]= 0.0; for (i=qh num_points; i--; ) { paraboloid= 0.0; for (k=0; k < qh hull_dim-1; k++) { paraboloid += *coord * *coord; infinity[k] += *coord; coord++; } *(coord++)= paraboloid; maximize_(maxboloid, paraboloid); } /* coord == infinity */ for (k=qh hull_dim-1; k--; ) *(coord++) /= qh num_points; *(coord++)= maxboloid * 1.1; qh num_points++; trace0((qh ferr, 9, "qh_projectinput: projected points to paraboloid for Delaunay\n")); }else if (qh DELAUNAY) /* !qh ATinfinity */ qh_setdelaunay(qh hull_dim, qh num_points, qh first_point); } /* projectinput */ /*--------------------------------- qh_projectpoints( project, n, points, numpoints, dim, newpoints, newdim ) project points/numpoints/dim to newpoints/newdim if project[k] == -1 delete dimension k if project[k] == 1 add dimension k by duplicating previous column n is size of project notes: newpoints may be points if only adding dimension at end design: check that 'project' and 'newdim' agree for each dimension if project == -1 skip dimension else determine start of column in newpoints determine start of column in points if project == +1, duplicate previous column copy dimension (column) from points to newpoints */ void qh_projectpoints(signed char *project, int n, realT *points, int numpoints, int dim, realT *newpoints, int newdim) { int testdim= dim, oldk=0, newk=0, i,j=0,k; realT *newp, *oldp; for (k=0; k < n; k++) testdim += project[k]; if (testdim != newdim) { qh_fprintf(qh ferr, 6018, "qhull internal error (qh_projectpoints): newdim %d should be %d after projection\n", newdim, testdim); qh_errexit(qh_ERRqhull, NULL, NULL); } for (j=0; j= dim) continue; oldp= points+oldk; }else oldp= points+oldk++; for (i=numpoints; i--; ) { *newp= *oldp; newp += newdim; oldp += dim; } } if (oldk >= dim) break; } trace1((qh ferr, 1004, "qh_projectpoints: projected %d points from dim %d to dim %d\n", numpoints, dim, newdim)); } /* projectpoints */ /*--------------------------------- qh_rotateinput( rows ) rotate input using row matrix input points given by qh first_point, num_points, hull_dim assumes rows[dim] is a scratch buffer if qh POINTSmalloc, overwrites input points, else mallocs a new array returns: rotated input sets qh POINTSmalloc design: see qh_rotatepoints */ void qh_rotateinput(realT **rows) { if (!qh POINTSmalloc) { qh first_point= qh_copypoints(qh first_point, qh num_points, qh hull_dim); qh POINTSmalloc= True; } qh_rotatepoints(qh first_point, qh num_points, qh hull_dim, rows); } /* rotateinput */ /*--------------------------------- qh_rotatepoints( points, numpoints, dim, row ) rotate numpoints points by a d-dim row matrix assumes rows[dim] is a scratch buffer returns: rotated points in place design: for each point for each coordinate use row[dim] to compute partial inner product for each coordinate rotate by partial inner product */ void qh_rotatepoints(realT *points, int numpoints, int dim, realT **row) { realT *point, *rowi, *coord= NULL, sum, *newval; int i,j,k; if (qh IStracing >= 1) qh_printmatrix(qh ferr, "qh_rotatepoints: rotate points by", row, dim, dim); for (point=points, j=numpoints; j--; point += dim) { newval= row[dim]; for (i=0; i < dim; i++) { rowi= row[i]; coord= point; for (sum=0.0, k=dim; k--; ) sum += *rowi++ * *coord++; *(newval++)= sum; } for (k=dim; k--; ) *(--coord)= *(--newval); } } /* rotatepoints */ /*--------------------------------- qh_scaleinput( ) scale input points using qh low_bound/high_bound input points given by qh first_point, num_points, hull_dim if qh POINTSmalloc, overwrites input points, else mallocs a new array returns: scales coordinates of points to low_bound[k], high_bound[k] sets qh POINTSmalloc design: see qh_scalepoints */ void qh_scaleinput(void) { if (!qh POINTSmalloc) { qh first_point= qh_copypoints(qh first_point, qh num_points, qh hull_dim); qh POINTSmalloc= True; } qh_scalepoints(qh first_point, qh num_points, qh hull_dim, qh lower_bound, qh upper_bound); } /* scaleinput */ /*--------------------------------- qh_scalelast( points, numpoints, dim, low, high, newhigh ) scale last coordinate to [0.0, newhigh], for Delaunay triangulation input points given by points, numpoints, dim returns: changes scale of last coordinate from [low, high] to [0.0, newhigh] overwrites last coordinate of each point saves low/high/newhigh in qh.last_low, etc. for qh_setdelaunay() notes: to reduce precision issues, qh_scalelast makes the last coordinate similar to other coordinates the last coordinate for Delaunay triangulation is the sum of squares of input coordinates note that the range [0.0, newwidth] is wrong for narrow distributions with large positive coordinates (e.g., [995933.64, 995963.48]) when called by qh_setdelaunay, low/high may not match the data passed to qh_setdelaunay design: compute scale and shift factors apply to last coordinate of each point */ void qh_scalelast(coordT *points, int numpoints, int dim, coordT low, coordT high, coordT newhigh) { realT scale, shift; coordT *coord, newlow; int i; boolT nearzero= False; newlow= 0.0; trace4((qh ferr, 4013, "qh_scalelast: scale last coordinate from [%2.2g, %2.2g] to [%2.2g, %2.2g]\n", low, high, newlow, newhigh)); qh last_low= low; qh last_high= high; qh last_newhigh= newhigh; scale= qh_divzero(newhigh - newlow, high - low, qh MINdenom_1, &nearzero); if (nearzero) { if (qh DELAUNAY) qh_fprintf(qh ferr, 6019, "qhull input error (qh_scalelast): can not scale last coordinate to [%4.4g, %4.4g]. Input is cocircular or cospherical. Use option 'Qz' to add a point at infinity.\n", newlow, newhigh); else qh_fprintf(qh ferr, 6020, "qhull input error (qh_scalelast): can not scale last coordinate to [%4.4g, %4.4g]. New bounds are too wide for compared to existing bounds [%4.4g, %4.4g] (width %4.4g)\n", newlow, newhigh, low, high, high-low); qh_errexit(qh_ERRinput, NULL, NULL); } shift= newlow - low * scale; coord= points + dim - 1; for (i=numpoints; i--; coord += dim) *coord= *coord * scale + shift; } /* scalelast */ /*--------------------------------- qh_scalepoints( points, numpoints, dim, newlows, newhighs ) scale points to new lowbound and highbound retains old bound when newlow= -REALmax or newhigh= +REALmax returns: scaled points overwrites old points design: for each coordinate compute current low and high bound compute scale and shift factors scale all points enforce new low and high bound for all points */ void qh_scalepoints(pointT *points, int numpoints, int dim, realT *newlows, realT *newhighs) { int i,k; realT shift, scale, *coord, low, high, newlow, newhigh, mincoord, maxcoord; boolT nearzero= False; for (k=0; k < dim; k++) { newhigh= newhighs[k]; newlow= newlows[k]; if (newhigh > REALmax/2 && newlow < -REALmax/2) continue; low= REALmax; high= -REALmax; for (i=numpoints, coord=points+k; i--; coord += dim) { minimize_(low, *coord); maximize_(high, *coord); } if (newhigh > REALmax/2) newhigh= high; if (newlow < -REALmax/2) newlow= low; if (qh DELAUNAY && k == dim-1 && newhigh < newlow) { qh_fprintf(qh ferr, 6021, "qhull input error: 'Qb%d' or 'QB%d' inverts paraboloid since high bound %.2g < low bound %.2g\n", k, k, newhigh, newlow); qh_errexit(qh_ERRinput, NULL, NULL); } scale= qh_divzero(newhigh - newlow, high - low, qh MINdenom_1, &nearzero); if (nearzero) { qh_fprintf(qh ferr, 6022, "qhull input error: %d'th dimension's new bounds [%2.2g, %2.2g] too wide for\nexisting bounds [%2.2g, %2.2g]\n", k, newlow, newhigh, low, high); qh_errexit(qh_ERRinput, NULL, NULL); } shift= (newlow * high - low * newhigh)/(high-low); coord= points+k; for (i=numpoints; i--; coord += dim) *coord= *coord * scale + shift; coord= points+k; if (newlow < newhigh) { mincoord= newlow; maxcoord= newhigh; }else { mincoord= newhigh; maxcoord= newlow; } for (i=numpoints; i--; coord += dim) { minimize_(*coord, maxcoord); /* because of roundoff error */ maximize_(*coord, mincoord); } trace0((qh ferr, 10, "qh_scalepoints: scaled %d'th coordinate [%2.2g, %2.2g] to [%.2g, %.2g] for %d points by %2.2g and shifted %2.2g\n", k, low, high, newlow, newhigh, numpoints, scale, shift)); } } /* scalepoints */ /*--------------------------------- qh_setdelaunay( dim, count, points ) project count points to dim-d paraboloid for Delaunay triangulation dim is one more than the dimension of the input set assumes dim is at least 3 (i.e., at least a 2-d Delaunay triangulation) points is a dim*count realT array. The first dim-1 coordinates are the coordinates of the first input point. array[dim] is the first coordinate of the second input point. array[2*dim] is the first coordinate of the third input point. if qh.last_low defined (i.e., 'Qbb' called qh_scalelast) calls qh_scalelast to scale the last coordinate the same as the other points returns: for each point sets point[dim-1] to sum of squares of coordinates scale points to 'Qbb' if needed notes: to project one point, use qh_setdelaunay(qh hull_dim, 1, point) Do not use options 'Qbk', 'QBk', or 'QbB' since they scale the coordinates after the original projection. */ void qh_setdelaunay(int dim, int count, pointT *points) { int i, k; coordT *coordp, coord; realT paraboloid; trace0((qh ferr, 11, "qh_setdelaunay: project %d points to paraboloid for Delaunay triangulation\n", count)); coordp= points; for (i=0; i < count; i++) { coord= *coordp++; paraboloid= coord*coord; for (k=dim-2; k--; ) { coord= *coordp++; paraboloid += coord*coord; } *coordp++= paraboloid; } if (qh last_low < REALmax/2) qh_scalelast(points, count, dim, qh last_low, qh last_high, qh last_newhigh); } /* setdelaunay */ /*--------------------------------- qh_sethalfspace( dim, coords, nextp, normal, offset, feasible ) set point to dual of halfspace relative to feasible point halfspace is normal coefficients and offset. returns: false and prints error if feasible point is outside of hull overwrites coordinates for point at dim coords nextp= next point (coords) does not call qh_errexit design: compute distance from feasible point to halfspace divide each normal coefficient by -dist */ boolT qh_sethalfspace(int dim, coordT *coords, coordT **nextp, coordT *normal, coordT *offset, coordT *feasible) { coordT *normp= normal, *feasiblep= feasible, *coordp= coords; realT dist; realT r; /*bug fix*/ int k; boolT zerodiv; dist= *offset; for (k=dim; k--; ) dist += *(normp++) * *(feasiblep++); if (dist > 0) goto LABELerroroutside; normp= normal; if (dist < -qh MINdenom) { for (k=dim; k--; ) *(coordp++)= *(normp++) / -dist; }else { for (k=dim; k--; ) { *(coordp++)= qh_divzero(*(normp++), -dist, qh MINdenom_1, &zerodiv); if (zerodiv) goto LABELerroroutside; } } *nextp= coordp; #ifndef qh_NOtrace if (qh IStracing >= 4) { qh_fprintf(qh ferr, 8021, "qh_sethalfspace: halfspace at offset %6.2g to point: ", *offset); for (k=dim, coordp=coords; k--; ) { r= *coordp++; qh_fprintf(qh ferr, 8022, " %6.2g", r); } qh_fprintf(qh ferr, 8023, "\n"); } #endif return True; LABELerroroutside: feasiblep= feasible; normp= normal; qh_fprintf(qh ferr, 6023, "qhull input error: feasible point is not clearly inside halfspace\nfeasible point: "); for (k=dim; k--; ) qh_fprintf(qh ferr, 8024, qh_REAL_1, r=*(feasiblep++)); qh_fprintf(qh ferr, 8025, "\n halfspace: "); for (k=dim; k--; ) qh_fprintf(qh ferr, 8026, qh_REAL_1, r=*(normp++)); qh_fprintf(qh ferr, 8027, "\n at offset: "); qh_fprintf(qh ferr, 8028, qh_REAL_1, *offset); qh_fprintf(qh ferr, 8029, " and distance: "); qh_fprintf(qh ferr, 8030, qh_REAL_1, dist); qh_fprintf(qh ferr, 8031, "\n"); return False; } /* sethalfspace */ /*--------------------------------- qh_sethalfspace_all( dim, count, halfspaces, feasible ) generate dual for halfspace intersection with feasible point array of count halfspaces each halfspace is normal coefficients followed by offset the origin is inside the halfspace if the offset is negative feasible is a point inside all halfspaces (http://www.qhull.org/html/qhalf.htm#notes) returns: malloc'd array of count X dim-1 points notes: call before qh_init_B or qh_initqhull_globals free memory when done unused/untested code: please email bradb@shore.net if this works ok for you if using option 'Fp', qh.feasible_point must be set (e.g., to 'feasible') qh feasible_point is a malloc'd array that is freed by qh_freebuffers. design: see qh_sethalfspace */ coordT *qh_sethalfspace_all(int dim, int count, coordT *halfspaces, pointT *feasible) { int i, newdim; pointT *newpoints; coordT *coordp, *normalp, *offsetp; trace0((qh ferr, 12, "qh_sethalfspace_all: compute dual for halfspace intersection\n")); newdim= dim - 1; if (!(newpoints= (coordT *)qh_malloc((size_t)(count * newdim) * sizeof(coordT)))){ qh_fprintf(qh ferr, 6024, "qhull error: insufficient memory to compute dual of %d halfspaces\n", count); qh_errexit(qh_ERRmem, NULL, NULL); } coordp= newpoints; normalp= halfspaces; for (i=0; i < count; i++) { offsetp= normalp + newdim; if (!qh_sethalfspace(newdim, coordp, &coordp, normalp, offsetp, feasible)) { qh_free(newpoints); /* feasible is not inside halfspace as reported by qh_sethalfspace */ qh_fprintf(qh ferr, 8032, "The halfspace was at index %d\n", i); qh_errexit(qh_ERRinput, NULL, NULL); } normalp= offsetp + 1; } return newpoints; } /* sethalfspace_all */ /*--------------------------------- qh_sharpnewfacets( ) returns: true if could be an acute angle (facets in different quadrants) notes: for qh_findbest design: for all facets on qh.newfacet_list if two facets are in different quadrants set issharp */ boolT qh_sharpnewfacets(void) { facetT *facet; boolT issharp= False; int *quadrant, k; quadrant= (int *)qh_memalloc(qh hull_dim * (int)sizeof(int)); FORALLfacet_(qh newfacet_list) { if (facet == qh newfacet_list) { for (k=qh hull_dim; k--; ) quadrant[ k]= (facet->normal[ k] > 0); }else { for (k=qh hull_dim; k--; ) { if (quadrant[ k] != (facet->normal[ k] > 0)) { issharp= True; break; } } } if (issharp) break; } qh_memfree(quadrant, qh hull_dim * (int)sizeof(int)); trace3((qh ferr, 3001, "qh_sharpnewfacets: %d\n", issharp)); return issharp; } /* sharpnewfacets */ /*--------------------------------- qh_vertex_bestdist( vertices ) qh_vertex_bestdist2( vertices, vertexp, vertexp2 ) return nearest distance between vertices optionally returns vertex and vertex2 notes: called by qh_partitioncoplanar, qh_mergefacet, qh_check_maxout, qh_checkpoint */ coordT qh_vertex_bestdist(setT *vertices) { vertexT *vertex, *vertex2; return qh_vertex_bestdist2(vertices, &vertex, &vertex2); } /* vertex_bestdist */ coordT qh_vertex_bestdist2(setT *vertices, vertexT **vertexp/*= NULL*/, vertexT **vertexp2/*= NULL*/) { vertexT *vertex, *vertexA, *bestvertex= NULL, *bestvertex2= NULL; coordT dist, bestdist= REALmax; int k, vertex_i, vertex_n; FOREACHvertex_i_(vertices) { for (k= vertex_i+1; k < vertex_n; k++) { vertexA= SETelemt_(vertices, k, vertexT); dist= qh_pointdist(vertex->point, vertexA->point, -qh hull_dim); if (dist < bestdist) { bestdist= dist; bestvertex= vertex; bestvertex2= vertexA; } } } *vertexp= bestvertex; *vertexp2= bestvertex2; return sqrt(bestdist); } /* vertex_bestdist */ /*--------------------------------- qh_voronoi_center( dim, points ) return Voronoi center for a set of points dim is the orginal dimension of the points gh.gm_matrix/qh.gm_row are scratch buffers returns: center as a temporary point (qh_memalloc) if non-simplicial, returns center for max simplex of points notes: only called by qh_facetcenter from Bowyer & Woodwark, A Programmer's Geometry, 1983, p. 65 design: if non-simplicial determine max simplex for points translate point0 of simplex to origin compute sum of squares of diagonal compute determinate compute Voronoi center (see Bowyer & Woodwark) */ pointT *qh_voronoi_center(int dim, setT *points) { pointT *point, **pointp, *point0; pointT *center= (pointT *)qh_memalloc(qh center_size); setT *simplex; int i, j, k, size= qh_setsize(points); coordT *gmcoord; realT *diffp, sum2, *sum2row, *sum2p, det, factor; boolT nearzero, infinite; if (size == dim+1) simplex= points; else if (size < dim+1) { qh_memfree(center, qh center_size); qh_fprintf(qh ferr, 6025, "qhull internal error (qh_voronoi_center): need at least %d points to construct a Voronoi center\n", dim+1); qh_errexit(qh_ERRqhull, NULL, NULL); simplex= points; /* never executed -- avoids warning */ }else { simplex= qh_settemp(dim+1); qh_maxsimplex(dim, points, NULL, 0, &simplex); } point0= SETfirstt_(simplex, pointT); gmcoord= qh gm_matrix; for (k=0; k < dim; k++) { qh gm_row[k]= gmcoord; FOREACHpoint_(simplex) { if (point != point0) *(gmcoord++)= point[k] - point0[k]; } } sum2row= gmcoord; for (i=0; i < dim; i++) { sum2= 0.0; for (k=0; k < dim; k++) { diffp= qh gm_row[k] + i; sum2 += *diffp * *diffp; } *(gmcoord++)= sum2; } det= qh_determinant(qh gm_row, dim, &nearzero); factor= qh_divzero(0.5, det, qh MINdenom, &infinite); if (infinite) { for (k=dim; k--; ) center[k]= qh_INFINITE; if (qh IStracing) qh_printpoints(qh ferr, "qh_voronoi_center: at infinity for ", simplex); }else { for (i=0; i < dim; i++) { gmcoord= qh gm_matrix; sum2p= sum2row; for (k=0; k < dim; k++) { qh gm_row[k]= gmcoord; if (k == i) { for (j=dim; j--; ) *(gmcoord++)= *sum2p++; }else { FOREACHpoint_(simplex) { if (point != point0) *(gmcoord++)= point[k] - point0[k]; } } } center[i]= qh_determinant(qh gm_row, dim, &nearzero)*factor + point0[i]; } #ifndef qh_NOtrace if (qh IStracing >= 3) { qh_fprintf(qh ferr, 3061, "qh_voronoi_center: det %2.2g factor %2.2g ", det, factor); qh_printmatrix(qh ferr, "center:", ¢er, 1, dim); if (qh IStracing >= 5) { qh_printpoints(qh ferr, "points", simplex); FOREACHpoint_(simplex) qh_fprintf(qh ferr, 8034, "p%d dist %.2g, ", qh_pointid(point), qh_pointdist(point, center, dim)); qh_fprintf(qh ferr, 8035, "\n"); } } #endif } if (simplex != points) qh_settempfree(&simplex); return center; } /* voronoi_center */