geom2.c
Go to the documentation of this file.
00001 /*<html><pre>  -<a                             href="qh-geom.htm"
00002   >-------------------------------</a><a name="TOP">-</a>
00003 
00004 
00005    geom2.c
00006    infrequently used geometric routines of qhull
00007 
00008    see qh-geom.htm and geom.h
00009 
00010    Copyright (c) 1993-2011 The Geometry Center.
00011    $Id: //main/2011/qhull/src/libqhull/geom2.c#2 $$Change: 1342 $
00012    $DateTime: 2011/03/07 21:55:47 $$Author: bbarber $
00013 
00014    frequently used code goes into geom.c
00015 */
00016 
00017 #include "qhull_a.h"
00018 
00019 /*================== functions in alphabetic order ============*/
00020 
00021 /*-<a                             href="qh-geom.htm#TOC"
00022   >-------------------------------</a><a name="copypoints">-</a>
00023 
00024   qh_copypoints( points, numpoints, dimension)
00025     return qh_malloc'd copy of points
00026 */
00027 coordT *qh_copypoints(coordT *points, int numpoints, int dimension) {
00028   int size;
00029   coordT *newpoints;
00030 
00031   size= numpoints * dimension * (int)sizeof(coordT);
00032   if (!(newpoints=(coordT*)qh_malloc((size_t)size))) {
00033     qh_fprintf(qh ferr, 6004, "qhull error: insufficient memory to copy %d points\n",
00034         numpoints);
00035     qh_errexit(qh_ERRmem, NULL, NULL);
00036   }
00037   memcpy((char *)newpoints, (char *)points, (size_t)size);
00038   return newpoints;
00039 } /* copypoints */
00040 
00041 /*-<a                             href="qh-geom.htm#TOC"
00042   >-------------------------------</a><a name="crossproduct">-</a>
00043 
00044   qh_crossproduct( dim, vecA, vecB, vecC )
00045     crossproduct of 2 dim vectors
00046     C= A x B
00047 
00048   notes:
00049     from Glasner, Graphics Gems I, p. 639
00050     only defined for dim==3
00051 */
00052 void qh_crossproduct(int dim, realT vecA[3], realT vecB[3], realT vecC[3]){
00053 
00054   if (dim == 3) {
00055     vecC[0]=   det2_(vecA[1], vecA[2],
00056                      vecB[1], vecB[2]);
00057     vecC[1]= - det2_(vecA[0], vecA[2],
00058                      vecB[0], vecB[2]);
00059     vecC[2]=   det2_(vecA[0], vecA[1],
00060                      vecB[0], vecB[1]);
00061   }
00062 } /* vcross */
00063 
00064 /*-<a                             href="qh-geom.htm#TOC"
00065   >-------------------------------</a><a name="determinant">-</a>
00066 
00067   qh_determinant( rows, dim, nearzero )
00068     compute signed determinant of a square matrix
00069     uses qh.NEARzero to test for degenerate matrices
00070 
00071   returns:
00072     determinant
00073     overwrites rows and the matrix
00074     if dim == 2 or 3
00075       nearzero iff determinant < qh NEARzero[dim-1]
00076       (!quite correct, not critical)
00077     if dim >= 4
00078       nearzero iff diagonal[k] < qh NEARzero[k]
00079 */
00080 realT qh_determinant(realT **rows, int dim, boolT *nearzero) {
00081   realT det=0;
00082   int i;
00083   boolT sign= False;
00084 
00085   *nearzero= False;
00086   if (dim < 2) {
00087     qh_fprintf(qh ferr, 6005, "qhull internal error (qh_determinate): only implemented for dimension >= 2\n");
00088     qh_errexit(qh_ERRqhull, NULL, NULL);
00089   }else if (dim == 2) {
00090     det= det2_(rows[0][0], rows[0][1],
00091                  rows[1][0], rows[1][1]);
00092     if (fabs_(det) < qh NEARzero[1])  /* not really correct, what should this be? */
00093       *nearzero= True;
00094   }else if (dim == 3) {
00095     det= det3_(rows[0][0], rows[0][1], rows[0][2],
00096                  rows[1][0], rows[1][1], rows[1][2],
00097                  rows[2][0], rows[2][1], rows[2][2]);
00098     if (fabs_(det) < qh NEARzero[2])  /* not really correct, what should this be? */
00099       *nearzero= True;
00100   }else {
00101     qh_gausselim(rows, dim, dim, &sign, nearzero);  /* if nearzero, diagonal still ok*/
00102     det= 1.0;
00103     for (i=dim; i--; )
00104       det *= (rows[i])[i];
00105     if (sign)
00106       det= -det;
00107   }
00108   return det;
00109 } /* determinant */
00110 
00111 /*-<a                             href="qh-geom.htm#TOC"
00112   >-------------------------------</a><a name="detjoggle">-</a>
00113 
00114   qh_detjoggle( points, numpoints, dimension )
00115     determine default max joggle for point array
00116       as qh_distround * qh_JOGGLEdefault
00117 
00118   returns:
00119     initial value for JOGGLEmax from points and REALepsilon
00120 
00121   notes:
00122     computes DISTround since qh_maxmin not called yet
00123     if qh SCALElast, last dimension will be scaled later to MAXwidth
00124 
00125     loop duplicated from qh_maxmin
00126 */
00127 realT qh_detjoggle(pointT *points, int numpoints, int dimension) {
00128   realT abscoord, distround, joggle, maxcoord, mincoord;
00129   pointT *point, *pointtemp;
00130   realT maxabs= -REALmax;
00131   realT sumabs= 0;
00132   realT maxwidth= 0;
00133   int k;
00134 
00135   for (k=0; k < dimension; k++) {
00136     if (qh SCALElast && k == dimension-1)
00137       abscoord= maxwidth;
00138     else if (qh DELAUNAY && k == dimension-1) /* will qh_setdelaunay() */
00139       abscoord= 2 * maxabs * maxabs;  /* may be low by qh hull_dim/2 */
00140     else {
00141       maxcoord= -REALmax;
00142       mincoord= REALmax;
00143       FORALLpoint_(points, numpoints) {
00144         maximize_(maxcoord, point[k]);
00145         minimize_(mincoord, point[k]);
00146       }
00147       maximize_(maxwidth, maxcoord-mincoord);
00148       abscoord= fmax_(maxcoord, -mincoord);
00149     }
00150     sumabs += abscoord;
00151     maximize_(maxabs, abscoord);
00152   } /* for k */
00153   distround= qh_distround(qh hull_dim, maxabs, sumabs);
00154   joggle= distround * qh_JOGGLEdefault;
00155   maximize_(joggle, REALepsilon * qh_JOGGLEdefault);
00156   trace2((qh ferr, 2001, "qh_detjoggle: joggle=%2.2g maxwidth=%2.2g\n", joggle, maxwidth));
00157   return joggle;
00158 } /* detjoggle */
00159 
00160 /*-<a                             href="qh-geom.htm#TOC"
00161   >-------------------------------</a><a name="detroundoff">-</a>
00162 
00163   qh_detroundoff()
00164     determine maximum roundoff errors from
00165       REALepsilon, REALmax, REALmin, qh.hull_dim, qh.MAXabs_coord,
00166       qh.MAXsumcoord, qh.MAXwidth, qh.MINdenom_1
00167 
00168     accounts for qh.SETroundoff, qh.RANDOMdist, qh MERGEexact
00169       qh.premerge_cos, qh.postmerge_cos, qh.premerge_centrum,
00170       qh.postmerge_centrum, qh.MINoutside,
00171       qh_RATIOnearinside, qh_COPLANARratio, qh_WIDEcoplanar
00172 
00173   returns:
00174     sets qh.DISTround, etc. (see below)
00175     appends precision constants to qh.qhull_options
00176 
00177   see:
00178     qh_maxmin() for qh.NEARzero
00179 
00180   design:
00181     determine qh.DISTround for distance computations
00182     determine minimum denominators for qh_divzero
00183     determine qh.ANGLEround for angle computations
00184     adjust qh.premerge_cos,... for roundoff error
00185     determine qh.ONEmerge for maximum error due to a single merge
00186     determine qh.NEARinside, qh.MAXcoplanar, qh.MINvisible,
00187       qh.MINoutside, qh.WIDEfacet
00188     initialize qh.max_vertex and qh.minvertex
00189 */
00190 void qh_detroundoff(void) {
00191 
00192   qh_option("_max-width", NULL, &qh MAXwidth);
00193   if (!qh SETroundoff) {
00194     qh DISTround= qh_distround(qh hull_dim, qh MAXabs_coord, qh MAXsumcoord);
00195     if (qh RANDOMdist)
00196       qh DISTround += qh RANDOMfactor * qh MAXabs_coord;
00197     qh_option("Error-roundoff", NULL, &qh DISTround);
00198   }
00199   qh MINdenom= qh MINdenom_1 * qh MAXabs_coord;
00200   qh MINdenom_1_2= sqrt(qh MINdenom_1 * qh hull_dim) ;  /* if will be normalized */
00201   qh MINdenom_2= qh MINdenom_1_2 * qh MAXabs_coord;
00202                                               /* for inner product */
00203   qh ANGLEround= 1.01 * qh hull_dim * REALepsilon;
00204   if (qh RANDOMdist)
00205     qh ANGLEround += qh RANDOMfactor;
00206   if (qh premerge_cos < REALmax/2) {
00207     qh premerge_cos -= qh ANGLEround;
00208     if (qh RANDOMdist)
00209       qh_option("Angle-premerge-with-random", NULL, &qh premerge_cos);
00210   }
00211   if (qh postmerge_cos < REALmax/2) {
00212     qh postmerge_cos -= qh ANGLEround;
00213     if (qh RANDOMdist)
00214       qh_option("Angle-postmerge-with-random", NULL, &qh postmerge_cos);
00215   }
00216   qh premerge_centrum += 2 * qh DISTround;    /*2 for centrum and distplane()*/
00217   qh postmerge_centrum += 2 * qh DISTround;
00218   if (qh RANDOMdist && (qh MERGEexact || qh PREmerge))
00219     qh_option("Centrum-premerge-with-random", NULL, &qh premerge_centrum);
00220   if (qh RANDOMdist && qh POSTmerge)
00221     qh_option("Centrum-postmerge-with-random", NULL, &qh postmerge_centrum);
00222   { /* compute ONEmerge, max vertex offset for merging simplicial facets */
00223     realT maxangle= 1.0, maxrho;
00224 
00225     minimize_(maxangle, qh premerge_cos);
00226     minimize_(maxangle, qh postmerge_cos);
00227     /* max diameter * sin theta + DISTround for vertex to its hyperplane */
00228     qh ONEmerge= sqrt((realT)qh hull_dim) * qh MAXwidth *
00229       sqrt(1.0 - maxangle * maxangle) + qh DISTround;
00230     maxrho= qh hull_dim * qh premerge_centrum + qh DISTround;
00231     maximize_(qh ONEmerge, maxrho);
00232     maxrho= qh hull_dim * qh postmerge_centrum + qh DISTround;
00233     maximize_(qh ONEmerge, maxrho);
00234     if (qh MERGING)
00235       qh_option("_one-merge", NULL, &qh ONEmerge);
00236   }
00237   qh NEARinside= qh ONEmerge * qh_RATIOnearinside; /* only used if qh KEEPnearinside */
00238   if (qh JOGGLEmax < REALmax/2 && (qh KEEPcoplanar || qh KEEPinside)) {
00239     realT maxdist;             /* adjust qh.NEARinside for joggle */
00240     qh KEEPnearinside= True;
00241     maxdist= sqrt((realT)qh hull_dim) * qh JOGGLEmax + qh DISTround;
00242     maxdist= 2*maxdist;        /* vertex and coplanar point can joggle in opposite directions */
00243     maximize_(qh NEARinside, maxdist);  /* must agree with qh_nearcoplanar() */
00244   }
00245   if (qh KEEPnearinside)
00246     qh_option("_near-inside", NULL, &qh NEARinside);
00247   if (qh JOGGLEmax < qh DISTround) {
00248     qh_fprintf(qh ferr, 6006, "qhull error: the joggle for 'QJn', %.2g, is below roundoff for distance computations, %.2g\n",
00249          qh JOGGLEmax, qh DISTround);
00250     qh_errexit(qh_ERRinput, NULL, NULL);
00251   }
00252   if (qh MINvisible > REALmax/2) {
00253     if (!qh MERGING)
00254       qh MINvisible= qh DISTround;
00255     else if (qh hull_dim <= 3)
00256       qh MINvisible= qh premerge_centrum;
00257     else
00258       qh MINvisible= qh_COPLANARratio * qh premerge_centrum;
00259     if (qh APPROXhull && qh MINvisible > qh MINoutside)
00260       qh MINvisible= qh MINoutside;
00261     qh_option("Visible-distance", NULL, &qh MINvisible);
00262   }
00263   if (qh MAXcoplanar > REALmax/2) {
00264     qh MAXcoplanar= qh MINvisible;
00265     qh_option("U-coplanar-distance", NULL, &qh MAXcoplanar);
00266   }
00267   if (!qh APPROXhull) {             /* user may specify qh MINoutside */
00268     qh MINoutside= 2 * qh MINvisible;
00269     if (qh premerge_cos < REALmax/2)
00270       maximize_(qh MINoutside, (1- qh premerge_cos) * qh MAXabs_coord);
00271     qh_option("Width-outside", NULL, &qh MINoutside);
00272   }
00273   qh WIDEfacet= qh MINoutside;
00274   maximize_(qh WIDEfacet, qh_WIDEcoplanar * qh MAXcoplanar);
00275   maximize_(qh WIDEfacet, qh_WIDEcoplanar * qh MINvisible);
00276   qh_option("_wide-facet", NULL, &qh WIDEfacet);
00277   if (qh MINvisible > qh MINoutside + 3 * REALepsilon
00278   && !qh BESToutside && !qh FORCEoutput)
00279     qh_fprintf(qh ferr, 7001, "qhull input warning: minimum visibility V%.2g is greater than \nminimum outside W%.2g.  Flipped facets are likely.\n",
00280              qh MINvisible, qh MINoutside);
00281   qh max_vertex= qh DISTround;
00282   qh min_vertex= -qh DISTround;
00283   /* numeric constants reported in printsummary */
00284 } /* detroundoff */
00285 
00286 /*-<a                             href="qh-geom.htm#TOC"
00287   >-------------------------------</a><a name="detsimplex">-</a>
00288 
00289   qh_detsimplex( apex, points, dim, nearzero )
00290     compute determinant of a simplex with point apex and base points
00291 
00292   returns:
00293      signed determinant and nearzero from qh_determinant
00294 
00295   notes:
00296      uses qh.gm_matrix/qh.gm_row (assumes they're big enough)
00297 
00298   design:
00299     construct qm_matrix by subtracting apex from points
00300     compute determinate
00301 */
00302 realT qh_detsimplex(pointT *apex, setT *points, int dim, boolT *nearzero) {
00303   pointT *coorda, *coordp, *gmcoord, *point, **pointp;
00304   coordT **rows;
00305   int k,  i=0;
00306   realT det;
00307 
00308   zinc_(Zdetsimplex);
00309   gmcoord= qh gm_matrix;
00310   rows= qh gm_row;
00311   FOREACHpoint_(points) {
00312     if (i == dim)
00313       break;
00314     rows[i++]= gmcoord;
00315     coordp= point;
00316     coorda= apex;
00317     for (k=dim; k--; )
00318       *(gmcoord++)= *coordp++ - *coorda++;
00319   }
00320   if (i < dim) {
00321     qh_fprintf(qh ferr, 6007, "qhull internal error (qh_detsimplex): #points %d < dimension %d\n",
00322                i, dim);
00323     qh_errexit(qh_ERRqhull, NULL, NULL);
00324   }
00325   det= qh_determinant(rows, dim, nearzero);
00326   trace2((qh ferr, 2002, "qh_detsimplex: det=%2.2g for point p%d, dim %d, nearzero? %d\n",
00327           det, qh_pointid(apex), dim, *nearzero));
00328   return det;
00329 } /* detsimplex */
00330 
00331 /*-<a                             href="qh-geom.htm#TOC"
00332   >-------------------------------</a><a name="distnorm">-</a>
00333 
00334   qh_distnorm( dim, point, normal, offset )
00335     return distance from point to hyperplane at normal/offset
00336 
00337   returns:
00338     dist
00339 
00340   notes:
00341     dist > 0 if point is outside of hyperplane
00342 
00343   see:
00344     qh_distplane in geom.c
00345 */
00346 realT qh_distnorm(int dim, pointT *point, pointT *normal, realT *offsetp) {
00347   coordT *normalp= normal, *coordp= point;
00348   realT dist;
00349   int k;
00350 
00351   dist= *offsetp;
00352   for (k=dim; k--; )
00353     dist += *(coordp++) * *(normalp++);
00354   return dist;
00355 } /* distnorm */
00356 
00357 /*-<a                             href="qh-geom.htm#TOC"
00358   >-------------------------------</a><a name="distround">-</a>
00359 
00360   qh_distround(dimension, maxabs, maxsumabs )
00361     compute maximum round-off error for a distance computation
00362       to a normalized hyperplane
00363     maxabs is the maximum absolute value of a coordinate
00364     maxsumabs is the maximum possible sum of absolute coordinate values
00365 
00366   returns:
00367     max dist round for REALepsilon
00368 
00369   notes:
00370     calculate roundoff error according to
00371     Lemma 3.2-1 of Golub and van Loan "Matrix Computation"
00372     use sqrt(dim) since one vector is normalized
00373       or use maxsumabs since one vector is < 1
00374 */
00375 realT qh_distround(int dimension, realT maxabs, realT maxsumabs) {
00376   realT maxdistsum, maxround;
00377 
00378   maxdistsum= sqrt((realT)dimension) * maxabs;
00379   minimize_( maxdistsum, maxsumabs);
00380   maxround= REALepsilon * (dimension * maxdistsum * 1.01 + maxabs);
00381               /* adds maxabs for offset */
00382   trace4((qh ferr, 4008, "qh_distround: %2.2g maxabs %2.2g maxsumabs %2.2g maxdistsum %2.2g\n",
00383                  maxround, maxabs, maxsumabs, maxdistsum));
00384   return maxround;
00385 } /* distround */
00386 
00387 /*-<a                             href="qh-geom.htm#TOC"
00388   >-------------------------------</a><a name="divzero">-</a>
00389 
00390   qh_divzero( numer, denom, mindenom1, zerodiv )
00391     divide by a number that's nearly zero
00392     mindenom1= minimum denominator for dividing into 1.0
00393 
00394   returns:
00395     quotient
00396     sets zerodiv and returns 0.0 if it would overflow
00397 
00398   design:
00399     if numer is nearly zero and abs(numer) < abs(denom)
00400       return numer/denom
00401     else if numer is nearly zero
00402       return 0 and zerodiv
00403     else if denom/numer non-zero
00404       return numer/denom
00405     else
00406       return 0 and zerodiv
00407 */
00408 realT qh_divzero(realT numer, realT denom, realT mindenom1, boolT *zerodiv) {
00409   realT temp, numerx, denomx;
00410 
00411 
00412   if (numer < mindenom1 && numer > -mindenom1) {
00413     numerx= fabs_(numer);
00414     denomx= fabs_(denom);
00415     if (numerx < denomx) {
00416       *zerodiv= False;
00417       return numer/denom;
00418     }else {
00419       *zerodiv= True;
00420       return 0.0;
00421     }
00422   }
00423   temp= denom/numer;
00424   if (temp > mindenom1 || temp < -mindenom1) {
00425     *zerodiv= False;
00426     return numer/denom;
00427   }else {
00428     *zerodiv= True;
00429     return 0.0;
00430   }
00431 } /* divzero */
00432 
00433 
00434 /*-<a                             href="qh-geom.htm#TOC"
00435   >-------------------------------</a><a name="facetarea">-</a>
00436 
00437   qh_facetarea( facet )
00438     return area for a facet
00439 
00440   notes:
00441     if non-simplicial,
00442       uses centrum to triangulate facet and sums the projected areas.
00443     if (qh DELAUNAY),
00444       computes projected area instead for last coordinate
00445     assumes facet->normal exists
00446     projecting tricoplanar facets to the hyperplane does not appear to make a difference
00447 
00448   design:
00449     if simplicial
00450       compute area
00451     else
00452       for each ridge
00453         compute area from centrum to ridge
00454     negate area if upper Delaunay facet
00455 */
00456 realT qh_facetarea(facetT *facet) {
00457   vertexT *apex;
00458   pointT *centrum;
00459   realT area= 0.0;
00460   ridgeT *ridge, **ridgep;
00461 
00462   if (facet->simplicial) {
00463     apex= SETfirstt_(facet->vertices, vertexT);
00464     area= qh_facetarea_simplex(qh hull_dim, apex->point, facet->vertices,
00465                     apex, facet->toporient, facet->normal, &facet->offset);
00466   }else {
00467     if (qh CENTERtype == qh_AScentrum)
00468       centrum= facet->center;
00469     else
00470       centrum= qh_getcentrum(facet);
00471     FOREACHridge_(facet->ridges)
00472       area += qh_facetarea_simplex(qh hull_dim, centrum, ridge->vertices,
00473                  NULL, (boolT)(ridge->top == facet),  facet->normal, &facet->offset);
00474     if (qh CENTERtype != qh_AScentrum)
00475       qh_memfree(centrum, qh normal_size);
00476   }
00477   if (facet->upperdelaunay && qh DELAUNAY)
00478     area= -area;  /* the normal should be [0,...,1] */
00479   trace4((qh ferr, 4009, "qh_facetarea: f%d area %2.2g\n", facet->id, area));
00480   return area;
00481 } /* facetarea */
00482 
00483 /*-<a                             href="qh-geom.htm#TOC"
00484   >-------------------------------</a><a name="facetarea_simplex">-</a>
00485 
00486   qh_facetarea_simplex( dim, apex, vertices, notvertex, toporient, normal, offset )
00487     return area for a simplex defined by
00488       an apex, a base of vertices, an orientation, and a unit normal
00489     if simplicial or tricoplanar facet,
00490       notvertex is defined and it is skipped in vertices
00491 
00492   returns:
00493     computes area of simplex projected to plane [normal,offset]
00494     returns 0 if vertex too far below plane (qh WIDEfacet)
00495       vertex can't be apex of tricoplanar facet
00496 
00497   notes:
00498     if (qh DELAUNAY),
00499       computes projected area instead for last coordinate
00500     uses qh gm_matrix/gm_row and qh hull_dim
00501     helper function for qh_facetarea
00502 
00503   design:
00504     if Notvertex
00505       translate simplex to apex
00506     else
00507       project simplex to normal/offset
00508       translate simplex to apex
00509     if Delaunay
00510       set last row/column to 0 with -1 on diagonal
00511     else
00512       set last row to Normal
00513     compute determinate
00514     scale and flip sign for area
00515 */
00516 realT qh_facetarea_simplex(int dim, coordT *apex, setT *vertices,
00517         vertexT *notvertex,  boolT toporient, coordT *normal, realT *offset) {
00518   pointT *coorda, *coordp, *gmcoord;
00519   coordT **rows, *normalp;
00520   int k,  i=0;
00521   realT area, dist;
00522   vertexT *vertex, **vertexp;
00523   boolT nearzero;
00524 
00525   gmcoord= qh gm_matrix;
00526   rows= qh gm_row;
00527   FOREACHvertex_(vertices) {
00528     if (vertex == notvertex)
00529       continue;
00530     rows[i++]= gmcoord;
00531     coorda= apex;
00532     coordp= vertex->point;
00533     normalp= normal;
00534     if (notvertex) {
00535       for (k=dim; k--; )
00536         *(gmcoord++)= *coordp++ - *coorda++;
00537     }else {
00538       dist= *offset;
00539       for (k=dim; k--; )
00540         dist += *coordp++ * *normalp++;
00541       if (dist < -qh WIDEfacet) {
00542         zinc_(Znoarea);
00543         return 0.0;
00544       }
00545       coordp= vertex->point;
00546       normalp= normal;
00547       for (k=dim; k--; )
00548         *(gmcoord++)= (*coordp++ - dist * *normalp++) - *coorda++;
00549     }
00550   }
00551   if (i != dim-1) {
00552     qh_fprintf(qh ferr, 6008, "qhull internal error (qh_facetarea_simplex): #points %d != dim %d -1\n",
00553                i, dim);
00554     qh_errexit(qh_ERRqhull, NULL, NULL);
00555   }
00556   rows[i]= gmcoord;
00557   if (qh DELAUNAY) {
00558     for (i=0; i < dim-1; i++)
00559       rows[i][dim-1]= 0.0;
00560     for (k=dim; k--; )
00561       *(gmcoord++)= 0.0;
00562     rows[dim-1][dim-1]= -1.0;
00563   }else {
00564     normalp= normal;
00565     for (k=dim; k--; )
00566       *(gmcoord++)= *normalp++;
00567   }
00568   zinc_(Zdetsimplex);
00569   area= qh_determinant(rows, dim, &nearzero);
00570   if (toporient)
00571     area= -area;
00572   area *= qh AREAfactor;
00573   trace4((qh ferr, 4010, "qh_facetarea_simplex: area=%2.2g for point p%d, toporient %d, nearzero? %d\n",
00574           area, qh_pointid(apex), toporient, nearzero));
00575   return area;
00576 } /* facetarea_simplex */
00577 
00578 /*-<a                             href="qh-geom.htm#TOC"
00579   >-------------------------------</a><a name="facetcenter">-</a>
00580 
00581   qh_facetcenter( vertices )
00582     return Voronoi center (Voronoi vertex) for a facet's vertices
00583 
00584   returns:
00585     return temporary point equal to the center
00586 
00587   see:
00588     qh_voronoi_center()
00589 */
00590 pointT *qh_facetcenter(setT *vertices) {
00591   setT *points= qh_settemp(qh_setsize(vertices));
00592   vertexT *vertex, **vertexp;
00593   pointT *center;
00594 
00595   FOREACHvertex_(vertices)
00596     qh_setappend(&points, vertex->point);
00597   center= qh_voronoi_center(qh hull_dim-1, points);
00598   qh_settempfree(&points);
00599   return center;
00600 } /* facetcenter */
00601 
00602 /*-<a                             href="qh-geom.htm#TOC"
00603   >-------------------------------</a><a name="findgooddist">-</a>
00604 
00605   qh_findgooddist( point, facetA, dist, facetlist )
00606     find best good facet visible for point from facetA
00607     assumes facetA is visible from point
00608 
00609   returns:
00610     best facet, i.e., good facet that is furthest from point
00611       distance to best facet
00612       NULL if none
00613 
00614     moves good, visible facets (and some other visible facets)
00615       to end of qh facet_list
00616 
00617   notes:
00618     uses qh visit_id
00619 
00620   design:
00621     initialize bestfacet if facetA is good
00622     move facetA to end of facetlist
00623     for each facet on facetlist
00624       for each unvisited neighbor of facet
00625         move visible neighbors to end of facetlist
00626         update best good neighbor
00627         if no good neighbors, update best facet
00628 */
00629 facetT *qh_findgooddist(pointT *point, facetT *facetA, realT *distp,
00630                facetT **facetlist) {
00631   realT bestdist= -REALmax, dist;
00632   facetT *neighbor, **neighborp, *bestfacet=NULL, *facet;
00633   boolT goodseen= False;
00634 
00635   if (facetA->good) {
00636     zzinc_(Zcheckpart);  /* calls from check_bestdist occur after print stats */
00637     qh_distplane(point, facetA, &bestdist);
00638     bestfacet= facetA;
00639     goodseen= True;
00640   }
00641   qh_removefacet(facetA);
00642   qh_appendfacet(facetA);
00643   *facetlist= facetA;
00644   facetA->visitid= ++qh visit_id;
00645   FORALLfacet_(*facetlist) {
00646     FOREACHneighbor_(facet) {
00647       if (neighbor->visitid == qh visit_id)
00648         continue;
00649       neighbor->visitid= qh visit_id;
00650       if (goodseen && !neighbor->good)
00651         continue;
00652       zzinc_(Zcheckpart);
00653       qh_distplane(point, neighbor, &dist);
00654       if (dist > 0) {
00655         qh_removefacet(neighbor);
00656         qh_appendfacet(neighbor);
00657         if (neighbor->good) {
00658           goodseen= True;
00659           if (dist > bestdist) {
00660             bestdist= dist;
00661             bestfacet= neighbor;
00662           }
00663         }
00664       }
00665     }
00666   }
00667   if (bestfacet) {
00668     *distp= bestdist;
00669     trace2((qh ferr, 2003, "qh_findgooddist: p%d is %2.2g above good facet f%d\n",
00670       qh_pointid(point), bestdist, bestfacet->id));
00671     return bestfacet;
00672   }
00673   trace4((qh ferr, 4011, "qh_findgooddist: no good facet for p%d above f%d\n",
00674       qh_pointid(point), facetA->id));
00675   return NULL;
00676 }  /* findgooddist */
00677 
00678 /*-<a                             href="qh-geom.htm#TOC"
00679   >-------------------------------</a><a name="getarea">-</a>
00680 
00681   qh_getarea( facetlist )
00682     set area of all facets in facetlist
00683     collect statistics
00684     nop if hasAreaVolume
00685 
00686   returns:
00687     sets qh totarea/totvol to total area and volume of convex hull
00688     for Delaunay triangulation, computes projected area of the lower or upper hull
00689       ignores upper hull if qh ATinfinity
00690 
00691   notes:
00692     could compute outer volume by expanding facet area by rays from interior
00693     the following attempt at perpendicular projection underestimated badly:
00694       qh.totoutvol += (-dist + facet->maxoutside + qh DISTround)
00695                             * area/ qh hull_dim;
00696   design:
00697     for each facet on facetlist
00698       compute facet->area
00699       update qh.totarea and qh.totvol
00700 */
00701 void qh_getarea(facetT *facetlist) {
00702   realT area;
00703   realT dist;
00704   facetT *facet;
00705 
00706   if (qh hasAreaVolume)
00707     return;
00708   if (qh REPORTfreq)
00709     qh_fprintf(qh ferr, 8020, "computing area of each facet and volume of the convex hull\n");
00710   else
00711     trace1((qh ferr, 1001, "qh_getarea: computing volume and area for each facet\n"));
00712   qh totarea= qh totvol= 0.0;
00713   FORALLfacet_(facetlist) {
00714     if (!facet->normal)
00715       continue;
00716     if (facet->upperdelaunay && qh ATinfinity)
00717       continue;
00718     if (!facet->isarea) {
00719       facet->f.area= qh_facetarea(facet);
00720       facet->isarea= True;
00721     }
00722     area= facet->f.area;
00723     if (qh DELAUNAY) {
00724       if (facet->upperdelaunay == qh UPPERdelaunay)
00725         qh totarea += area;
00726     }else {
00727       qh totarea += area;
00728       qh_distplane(qh interior_point, facet, &dist);
00729       qh totvol += -dist * area/ qh hull_dim;
00730     }
00731     if (qh PRINTstatistics) {
00732       wadd_(Wareatot, area);
00733       wmax_(Wareamax, area);
00734       wmin_(Wareamin, area);
00735     }
00736   }
00737   qh hasAreaVolume= True;
00738 } /* getarea */
00739 
00740 /*-<a                             href="qh-geom.htm#TOC"
00741   >-------------------------------</a><a name="gram_schmidt">-</a>
00742 
00743   qh_gram_schmidt( dim, row )
00744     implements Gram-Schmidt orthogonalization by rows
00745 
00746   returns:
00747     false if zero norm
00748     overwrites rows[dim][dim]
00749 
00750   notes:
00751     see Golub & van Loan Algorithm 6.2-2
00752     overflow due to small divisors not handled
00753 
00754   design:
00755     for each row
00756       compute norm for row
00757       if non-zero, normalize row
00758       for each remaining rowA
00759         compute inner product of row and rowA
00760         reduce rowA by row * inner product
00761 */
00762 boolT qh_gram_schmidt(int dim, realT **row) {
00763   realT *rowi, *rowj, norm;
00764   int i, j, k;
00765 
00766   for (i=0; i < dim; i++) {
00767     rowi= row[i];
00768     for (norm= 0.0, k= dim; k--; rowi++)
00769       norm += *rowi * *rowi;
00770     norm= sqrt(norm);
00771     wmin_(Wmindenom, norm);
00772     if (norm == 0.0)  /* either 0 or overflow due to sqrt */
00773       return False;
00774     for (k=dim; k--; )
00775       *(--rowi) /= norm;
00776     for (j=i+1; j < dim; j++) {
00777       rowj= row[j];
00778       for (norm= 0.0, k=dim; k--; )
00779         norm += *rowi++ * *rowj++;
00780       for (k=dim; k--; )
00781         *(--rowj) -= *(--rowi) * norm;
00782     }
00783   }
00784   return True;
00785 } /* gram_schmidt */
00786 
00787 
00788 /*-<a                             href="qh-geom.htm#TOC"
00789   >-------------------------------</a><a name="inthresholds">-</a>
00790 
00791   qh_inthresholds( normal, angle )
00792     return True if normal within qh.lower_/upper_threshold
00793 
00794   returns:
00795     estimate of angle by summing of threshold diffs
00796       angle may be NULL
00797       smaller "angle" is better
00798 
00799   notes:
00800     invalid if qh.SPLITthresholds
00801 
00802   see:
00803     qh.lower_threshold in qh_initbuild()
00804     qh_initthresholds()
00805 
00806   design:
00807     for each dimension
00808       test threshold
00809 */
00810 boolT qh_inthresholds(coordT *normal, realT *angle) {
00811   boolT within= True;
00812   int k;
00813   realT threshold;
00814 
00815   if (angle)
00816     *angle= 0.0;
00817   for (k=0; k < qh hull_dim; k++) {
00818     threshold= qh lower_threshold[k];
00819     if (threshold > -REALmax/2) {
00820       if (normal[k] < threshold)
00821         within= False;
00822       if (angle) {
00823         threshold -= normal[k];
00824         *angle += fabs_(threshold);
00825       }
00826     }
00827     if (qh upper_threshold[k] < REALmax/2) {
00828       threshold= qh upper_threshold[k];
00829       if (normal[k] > threshold)
00830         within= False;
00831       if (angle) {
00832         threshold -= normal[k];
00833         *angle += fabs_(threshold);
00834       }
00835     }
00836   }
00837   return within;
00838 } /* inthresholds */
00839 
00840 
00841 /*-<a                             href="qh-geom.htm#TOC"
00842   >-------------------------------</a><a name="joggleinput">-</a>
00843 
00844   qh_joggleinput()
00845     randomly joggle input to Qhull by qh.JOGGLEmax
00846     initial input is qh.first_point/qh.num_points of qh.hull_dim
00847       repeated calls use qh.input_points/qh.num_points
00848 
00849   returns:
00850     joggles points at qh.first_point/qh.num_points
00851     copies data to qh.input_points/qh.input_malloc if first time
00852     determines qh.JOGGLEmax if it was zero
00853     if qh.DELAUNAY
00854       computes the Delaunay projection of the joggled points
00855 
00856   notes:
00857     if qh.DELAUNAY, unnecessarily joggles the last coordinate
00858     the initial 'QJn' may be set larger than qh_JOGGLEmaxincrease
00859 
00860   design:
00861     if qh.DELAUNAY
00862       set qh.SCALElast for reduced precision errors
00863     if first call
00864       initialize qh.input_points to the original input points
00865       if qh.JOGGLEmax == 0
00866         determine default qh.JOGGLEmax
00867     else
00868       increase qh.JOGGLEmax according to qh.build_cnt
00869     joggle the input by adding a random number in [-qh.JOGGLEmax,qh.JOGGLEmax]
00870     if qh.DELAUNAY
00871       sets the Delaunay projection
00872 */
00873 void qh_joggleinput(void) {
00874   int i, seed, size;
00875   coordT *coordp, *inputp;
00876   realT randr, randa, randb;
00877 
00878   if (!qh input_points) { /* first call */
00879     qh input_points= qh first_point;
00880     qh input_malloc= qh POINTSmalloc;
00881     size= qh num_points * qh hull_dim * sizeof(coordT);
00882     if (!(qh first_point=(coordT*)qh_malloc((size_t)size))) {
00883       qh_fprintf(qh ferr, 6009, "qhull error: insufficient memory to joggle %d points\n",
00884           qh num_points);
00885       qh_errexit(qh_ERRmem, NULL, NULL);
00886     }
00887     qh POINTSmalloc= True;
00888     if (qh JOGGLEmax == 0.0) {
00889       qh JOGGLEmax= qh_detjoggle(qh input_points, qh num_points, qh hull_dim);
00890       qh_option("QJoggle", NULL, &qh JOGGLEmax);
00891     }
00892   }else {                 /* repeated call */
00893     if (!qh RERUN && qh build_cnt > qh_JOGGLEretry) {
00894       if (((qh build_cnt-qh_JOGGLEretry-1) % qh_JOGGLEagain) == 0) {
00895         realT maxjoggle= qh MAXwidth * qh_JOGGLEmaxincrease;
00896         if (qh JOGGLEmax < maxjoggle) {
00897           qh JOGGLEmax *= qh_JOGGLEincrease;
00898           minimize_(qh JOGGLEmax, maxjoggle);
00899         }
00900       }
00901     }
00902     qh_option("QJoggle", NULL, &qh JOGGLEmax);
00903   }
00904   if (qh build_cnt > 1 && qh JOGGLEmax > fmax_(qh MAXwidth/4, 0.1)) {
00905       qh_fprintf(qh ferr, 6010, "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",
00906                 qh JOGGLEmax);
00907       qh_errexit(qh_ERRqhull, NULL, NULL);
00908   }
00909   /* for some reason, using qh ROTATErandom and qh_RANDOMseed does not repeat the run. Use 'TRn' instead */
00910   seed= qh_RANDOMint;
00911   qh_option("_joggle-seed", &seed, NULL);
00912   trace0((qh ferr, 6, "qh_joggleinput: joggle input by %2.2g with seed %d\n",
00913     qh JOGGLEmax, seed));
00914   inputp= qh input_points;
00915   coordp= qh first_point;
00916   randa= 2.0 * qh JOGGLEmax/qh_RANDOMmax;
00917   randb= -qh JOGGLEmax;
00918   size= qh num_points * qh hull_dim;
00919   for (i=size; i--; ) {
00920     randr= qh_RANDOMint;
00921     *(coordp++)= *(inputp++) + (randr * randa + randb);
00922   }
00923   if (qh DELAUNAY) {
00924     qh last_low= qh last_high= qh last_newhigh= REALmax;
00925     qh_setdelaunay(qh hull_dim, qh num_points, qh first_point);
00926   }
00927 } /* joggleinput */
00928 
00929 /*-<a                             href="qh-geom.htm#TOC"
00930   >-------------------------------</a><a name="maxabsval">-</a>
00931 
00932   qh_maxabsval( normal, dim )
00933     return pointer to maximum absolute value of a dim vector
00934     returns NULL if dim=0
00935 */
00936 realT *qh_maxabsval(realT *normal, int dim) {
00937   realT maxval= -REALmax;
00938   realT *maxp= NULL, *colp, absval;
00939   int k;
00940 
00941   for (k=dim, colp= normal; k--; colp++) {
00942     absval= fabs_(*colp);
00943     if (absval > maxval) {
00944       maxval= absval;
00945       maxp= colp;
00946     }
00947   }
00948   return maxp;
00949 } /* maxabsval */
00950 
00951 
00952 /*-<a                             href="qh-geom.htm#TOC"
00953   >-------------------------------</a><a name="maxmin">-</a>
00954 
00955   qh_maxmin( points, numpoints, dimension )
00956     return max/min points for each dimension
00957     determine max and min coordinates
00958 
00959   returns:
00960     returns a temporary set of max and min points
00961       may include duplicate points. Does not include qh.GOODpoint
00962     sets qh.NEARzero, qh.MAXabs_coord, qh.MAXsumcoord, qh.MAXwidth
00963          qh.MAXlastcoord, qh.MINlastcoord
00964     initializes qh.max_outside, qh.min_vertex, qh.WAScoplanar, qh.ZEROall_ok
00965 
00966   notes:
00967     loop duplicated in qh_detjoggle()
00968 
00969   design:
00970     initialize global precision variables
00971     checks definition of REAL...
00972     for each dimension
00973       for each point
00974         collect maximum and minimum point
00975       collect maximum of maximums and minimum of minimums
00976       determine qh.NEARzero for Gaussian Elimination
00977 */
00978 setT *qh_maxmin(pointT *points, int numpoints, int dimension) {
00979   int k;
00980   realT maxcoord, temp;
00981   pointT *minimum, *maximum, *point, *pointtemp;
00982   setT *set;
00983 
00984   qh max_outside= 0.0;
00985   qh MAXabs_coord= 0.0;
00986   qh MAXwidth= -REALmax;
00987   qh MAXsumcoord= 0.0;
00988   qh min_vertex= 0.0;
00989   qh WAScoplanar= False;
00990   if (qh ZEROcentrum)
00991     qh ZEROall_ok= True;
00992   if (REALmin < REALepsilon && REALmin < REALmax && REALmin > -REALmax
00993   && REALmax > 0.0 && -REALmax < 0.0)
00994     ; /* all ok */
00995   else {
00996     qh_fprintf(qh ferr, 6011, "qhull error: floating point constants in user.h are wrong\n\
00997 REALepsilon %g REALmin %g REALmax %g -REALmax %g\n",
00998              REALepsilon, REALmin, REALmax, -REALmax);
00999     qh_errexit(qh_ERRinput, NULL, NULL);
01000   }
01001   set= qh_settemp(2*dimension);
01002   for (k=0; k < dimension; k++) {
01003     if (points == qh GOODpointp)
01004       minimum= maximum= points + dimension;
01005     else
01006       minimum= maximum= points;
01007     FORALLpoint_(points, numpoints) {
01008       if (point == qh GOODpointp)
01009         continue;
01010       if (maximum[k] < point[k])
01011         maximum= point;
01012       else if (minimum[k] > point[k])
01013         minimum= point;
01014     }
01015     if (k == dimension-1) {
01016       qh MINlastcoord= minimum[k];
01017       qh MAXlastcoord= maximum[k];
01018     }
01019     if (qh SCALElast && k == dimension-1)
01020       maxcoord= qh MAXwidth;
01021     else {
01022       maxcoord= fmax_(maximum[k], -minimum[k]);
01023       if (qh GOODpointp) {
01024         temp= fmax_(qh GOODpointp[k], -qh GOODpointp[k]);
01025         maximize_(maxcoord, temp);
01026       }
01027       temp= maximum[k] - minimum[k];
01028       maximize_(qh MAXwidth, temp);
01029     }
01030     maximize_(qh MAXabs_coord, maxcoord);
01031     qh MAXsumcoord += maxcoord;
01032     qh_setappend(&set, maximum);
01033     qh_setappend(&set, minimum);
01034     /* calculation of qh NEARzero is based on error formula 4.4-13 of
01035        Golub & van Loan, authors say n^3 can be ignored and 10 be used in
01036        place of rho */
01037     qh NEARzero[k]= 80 * qh MAXsumcoord * REALepsilon;
01038   }
01039   if (qh IStracing >=1)
01040     qh_printpoints(qh ferr, "qh_maxmin: found the max and min points(by dim):", set);
01041   return(set);
01042 } /* maxmin */
01043 
01044 /*-<a                             href="qh-geom.htm#TOC"
01045   >-------------------------------</a><a name="maxouter">-</a>
01046 
01047   qh_maxouter()
01048     return maximum distance from facet to outer plane
01049     normally this is qh.max_outside+qh.DISTround
01050     does not include qh.JOGGLEmax
01051 
01052   see:
01053     qh_outerinner()
01054 
01055   notes:
01056     need to add another qh.DISTround if testing actual point with computation
01057 
01058   for joggle:
01059     qh_setfacetplane() updated qh.max_outer for Wnewvertexmax (max distance to vertex)
01060     need to use Wnewvertexmax since could have a coplanar point for a high
01061       facet that is replaced by a low facet
01062     need to add qh.JOGGLEmax if testing input points
01063 */
01064 realT qh_maxouter(void) {
01065   realT dist;
01066 
01067   dist= fmax_(qh max_outside, qh DISTround);
01068   dist += qh DISTround;
01069   trace4((qh ferr, 4012, "qh_maxouter: max distance from facet to outer plane is %2.2g max_outside is %2.2g\n", dist, qh max_outside));
01070   return dist;
01071 } /* maxouter */
01072 
01073 /*-<a                             href="qh-geom.htm#TOC"
01074   >-------------------------------</a><a name="maxsimplex">-</a>
01075 
01076   qh_maxsimplex( dim, maxpoints, points, numpoints, simplex )
01077     determines maximum simplex for a set of points
01078     starts from points already in simplex
01079     skips qh.GOODpointp (assumes that it isn't in maxpoints)
01080 
01081   returns:
01082     simplex with dim+1 points
01083 
01084   notes:
01085     assumes at least pointsneeded points in points
01086     maximizes determinate for x,y,z,w, etc.
01087     uses maxpoints as long as determinate is clearly non-zero
01088 
01089   design:
01090     initialize simplex with at least two points
01091       (find points with max or min x coordinate)
01092     for each remaining dimension
01093       add point that maximizes the determinate
01094         (use points from maxpoints first)
01095 */
01096 void qh_maxsimplex(int dim, setT *maxpoints, pointT *points, int numpoints, setT **simplex) {
01097   pointT *point, **pointp, *pointtemp, *maxpoint, *minx=NULL, *maxx=NULL;
01098   boolT nearzero, maxnearzero= False;
01099   int k, sizinit;
01100   realT maxdet= -REALmax, det, mincoord= REALmax, maxcoord= -REALmax;
01101 
01102   sizinit= qh_setsize(*simplex);
01103   if (sizinit < 2) {
01104     if (qh_setsize(maxpoints) >= 2) {
01105       FOREACHpoint_(maxpoints) {
01106         if (maxcoord < point[0]) {
01107           maxcoord= point[0];
01108           maxx= point;
01109         }
01110         if (mincoord > point[0]) {
01111           mincoord= point[0];
01112           minx= point;
01113         }
01114       }
01115     }else {
01116       FORALLpoint_(points, numpoints) {
01117         if (point == qh GOODpointp)
01118           continue;
01119         if (maxcoord < point[0]) {
01120           maxcoord= point[0];
01121           maxx= point;
01122         }
01123         if (mincoord > point[0]) {
01124           mincoord= point[0];
01125           minx= point;
01126         }
01127       }
01128     }
01129     qh_setunique(simplex, minx);
01130     if (qh_setsize(*simplex) < 2)
01131       qh_setunique(simplex, maxx);
01132     sizinit= qh_setsize(*simplex);
01133     if (sizinit < 2) {
01134       qh_precision("input has same x coordinate");
01135       if (zzval_(Zsetplane) > qh hull_dim+1) {
01136         qh_fprintf(qh ferr, 6012, "qhull precision error (qh_maxsimplex for voronoi_center):\n%d points with the same x coordinate.\n",
01137                  qh_setsize(maxpoints)+numpoints);
01138         qh_errexit(qh_ERRprec, NULL, NULL);
01139       }else {
01140         qh_fprintf(qh ferr, 6013, "qhull input error: input is less than %d-dimensional since it has the same x coordinate\n", qh hull_dim);
01141         qh_errexit(qh_ERRinput, NULL, NULL);
01142       }
01143     }
01144   }
01145   for (k=sizinit; k < dim+1; k++) {
01146     maxpoint= NULL;
01147     maxdet= -REALmax;
01148     FOREACHpoint_(maxpoints) {
01149       if (!qh_setin(*simplex, point)) {
01150         det= qh_detsimplex(point, *simplex, k, &nearzero);
01151         if ((det= fabs_(det)) > maxdet) {
01152           maxdet= det;
01153           maxpoint= point;
01154           maxnearzero= nearzero;
01155         }
01156       }
01157     }
01158     if (!maxpoint || maxnearzero) {
01159       zinc_(Zsearchpoints);
01160       if (!maxpoint) {
01161         trace0((qh ferr, 7, "qh_maxsimplex: searching all points for %d-th initial vertex.\n", k+1));
01162       }else {
01163         trace0((qh ferr, 8, "qh_maxsimplex: searching all points for %d-th initial vertex, better than p%d det %2.2g\n",
01164                 k+1, qh_pointid(maxpoint), maxdet));
01165       }
01166       FORALLpoint_(points, numpoints) {
01167         if (point == qh GOODpointp)
01168           continue;
01169         if (!qh_setin(*simplex, point)) {
01170           det= qh_detsimplex(point, *simplex, k, &nearzero);
01171           if ((det= fabs_(det)) > maxdet) {
01172             maxdet= det;
01173             maxpoint= point;
01174             maxnearzero= nearzero;
01175           }
01176         }
01177       }
01178     } /* !maxpoint */
01179     if (!maxpoint) {
01180       qh_fprintf(qh ferr, 6014, "qhull internal error (qh_maxsimplex): not enough points available\n");
01181       qh_errexit(qh_ERRqhull, NULL, NULL);
01182     }
01183     qh_setappend(simplex, maxpoint);
01184     trace1((qh ferr, 1002, "qh_maxsimplex: selected point p%d for %d`th initial vertex, det=%2.2g\n",
01185             qh_pointid(maxpoint), k+1, maxdet));
01186   } /* k */
01187 } /* maxsimplex */
01188 
01189 /*-<a                             href="qh-geom.htm#TOC"
01190   >-------------------------------</a><a name="minabsval">-</a>
01191 
01192   qh_minabsval( normal, dim )
01193     return minimum absolute value of a dim vector
01194 */
01195 realT qh_minabsval(realT *normal, int dim) {
01196   realT minval= 0;
01197   realT maxval= 0;
01198   realT *colp;
01199   int k;
01200 
01201   for (k=dim, colp=normal; k--; colp++) {
01202     maximize_(maxval, *colp);
01203     minimize_(minval, *colp);
01204   }
01205   return fmax_(maxval, -minval);
01206 } /* minabsval */
01207 
01208 
01209 /*-<a                             href="qh-geom.htm#TOC"
01210   >-------------------------------</a><a name="mindiff">-</a>
01211 
01212   qh_mindif ( vecA, vecB, dim )
01213     return index of min abs. difference of two vectors
01214 */
01215 int qh_mindiff(realT *vecA, realT *vecB, int dim) {
01216   realT mindiff= REALmax, diff;
01217   realT *vecAp= vecA, *vecBp= vecB;
01218   int k, mink= 0;
01219 
01220   for (k=0; k < dim; k++) {
01221     diff= *vecAp++ - *vecBp++;
01222     diff= fabs_(diff);
01223     if (diff < mindiff) {
01224       mindiff= diff;
01225       mink= k;
01226     }
01227   }
01228   return mink;
01229 } /* mindiff */
01230 
01231 
01232 
01233 /*-<a                             href="qh-geom.htm#TOC"
01234   >-------------------------------</a><a name="orientoutside">-</a>
01235 
01236   qh_orientoutside( facet  )
01237     make facet outside oriented via qh.interior_point
01238 
01239   returns:
01240     True if facet reversed orientation.
01241 */
01242 boolT qh_orientoutside(facetT *facet) {
01243   int k;
01244   realT dist;
01245 
01246   qh_distplane(qh interior_point, facet, &dist);
01247   if (dist > 0) {
01248     for (k=qh hull_dim; k--; )
01249       facet->normal[k]= -facet->normal[k];
01250     facet->offset= -facet->offset;
01251     return True;
01252   }
01253   return False;
01254 } /* orientoutside */
01255 
01256 /*-<a                             href="qh-geom.htm#TOC"
01257   >-------------------------------</a><a name="outerinner">-</a>
01258 
01259   qh_outerinner( facet, outerplane, innerplane  )
01260     if facet and qh.maxoutdone (i.e., qh_check_maxout)
01261       returns outer and inner plane for facet
01262     else
01263       returns maximum outer and inner plane
01264     accounts for qh.JOGGLEmax
01265 
01266   see:
01267     qh_maxouter(), qh_check_bestdist(), qh_check_points()
01268 
01269   notes:
01270     outerplaner or innerplane may be NULL
01271     facet is const
01272     Does not error (QhullFacet)
01273 
01274     includes qh.DISTround for actual points
01275     adds another qh.DISTround if testing with floating point arithmetic
01276 */
01277 void qh_outerinner(facetT *facet, realT *outerplane, realT *innerplane) {
01278   realT dist, mindist;
01279   vertexT *vertex, **vertexp;
01280 
01281   if (outerplane) {
01282     if (!qh_MAXoutside || !facet || !qh maxoutdone) {
01283       *outerplane= qh_maxouter();       /* includes qh.DISTround */
01284     }else { /* qh_MAXoutside ... */
01285 #if qh_MAXoutside
01286       *outerplane= facet->maxoutside + qh DISTround;
01287 #endif
01288 
01289     }
01290     if (qh JOGGLEmax < REALmax/2)
01291       *outerplane += qh JOGGLEmax * sqrt((realT)qh hull_dim);
01292   }
01293   if (innerplane) {
01294     if (facet) {
01295       mindist= REALmax;
01296       FOREACHvertex_(facet->vertices) {
01297         zinc_(Zdistio);
01298         qh_distplane(vertex->point, facet, &dist);
01299         minimize_(mindist, dist);
01300       }
01301       *innerplane= mindist - qh DISTround;
01302     }else
01303       *innerplane= qh min_vertex - qh DISTround;
01304     if (qh JOGGLEmax < REALmax/2)
01305       *innerplane -= qh JOGGLEmax * sqrt((realT)qh hull_dim);
01306   }
01307 } /* outerinner */
01308 
01309 /*-<a                             href="qh-geom.htm#TOC"
01310   >-------------------------------</a><a name="pointdist">-</a>
01311 
01312   qh_pointdist( point1, point2, dim )
01313     return distance between two points
01314 
01315   notes:
01316     returns distance squared if 'dim' is negative
01317 */
01318 coordT qh_pointdist(pointT *point1, pointT *point2, int dim) {
01319   coordT dist, diff;
01320   int k;
01321 
01322   dist= 0.0;
01323   for (k= (dim > 0 ? dim : -dim); k--; ) {
01324     diff= *point1++ - *point2++;
01325     dist += diff * diff;
01326   }
01327   if (dim > 0)
01328     return(sqrt(dist));
01329   return dist;
01330 } /* pointdist */
01331 
01332 
01333 /*-<a                             href="qh-geom.htm#TOC"
01334   >-------------------------------</a><a name="printmatrix">-</a>
01335 
01336   qh_printmatrix( fp, string, rows, numrow, numcol )
01337     print matrix to fp given by row vectors
01338     print string as header
01339 
01340   notes:
01341     print a vector by qh_printmatrix(fp, "", &vect, 1, len)
01342 */
01343 void qh_printmatrix(FILE *fp, const char *string, realT **rows, int numrow, int numcol) {
01344   realT *rowp;
01345   realT r; /*bug fix*/
01346   int i,k;
01347 
01348   qh_fprintf(fp, 9001, "%s\n", string);
01349   for (i=0; i < numrow; i++) {
01350     rowp= rows[i];
01351     for (k=0; k < numcol; k++) {
01352       r= *rowp++;
01353       qh_fprintf(fp, 9002, "%6.3g ", r);
01354     }
01355     qh_fprintf(fp, 9003, "\n");
01356   }
01357 } /* printmatrix */
01358 
01359 
01360 /*-<a                             href="qh-geom.htm#TOC"
01361   >-------------------------------</a><a name="printpoints">-</a>
01362 
01363   qh_printpoints( fp, string, points )
01364     print pointids to fp for a set of points
01365     if string, prints string and 'p' point ids
01366 */
01367 void qh_printpoints(FILE *fp, const char *string, setT *points) {
01368   pointT *point, **pointp;
01369 
01370   if (string) {
01371     qh_fprintf(fp, 9004, "%s", string);
01372     FOREACHpoint_(points)
01373       qh_fprintf(fp, 9005, " p%d", qh_pointid(point));
01374     qh_fprintf(fp, 9006, "\n");
01375   }else {
01376     FOREACHpoint_(points)
01377       qh_fprintf(fp, 9007, " %d", qh_pointid(point));
01378     qh_fprintf(fp, 9008, "\n");
01379   }
01380 } /* printpoints */
01381 
01382 
01383 /*-<a                             href="qh-geom.htm#TOC"
01384   >-------------------------------</a><a name="projectinput">-</a>
01385 
01386   qh_projectinput()
01387     project input points using qh.lower_bound/upper_bound and qh DELAUNAY
01388     if qh.lower_bound[k]=qh.upper_bound[k]= 0,
01389       removes dimension k
01390     if halfspace intersection
01391       removes dimension k from qh.feasible_point
01392     input points in qh first_point, num_points, input_dim
01393 
01394   returns:
01395     new point array in qh first_point of qh hull_dim coordinates
01396     sets qh POINTSmalloc
01397     if qh DELAUNAY
01398       projects points to paraboloid
01399       lowbound/highbound is also projected
01400     if qh ATinfinity
01401       adds point "at-infinity"
01402     if qh POINTSmalloc
01403       frees old point array
01404 
01405   notes:
01406     checks that qh.hull_dim agrees with qh.input_dim, PROJECTinput, and DELAUNAY
01407 
01408 
01409   design:
01410     sets project[k] to -1 (delete), 0 (keep), 1 (add for Delaunay)
01411     determines newdim and newnum for qh hull_dim and qh num_points
01412     projects points to newpoints
01413     projects qh.lower_bound to itself
01414     projects qh.upper_bound to itself
01415     if qh DELAUNAY
01416       if qh ATINFINITY
01417         projects points to paraboloid
01418         computes "infinity" point as vertex average and 10% above all points
01419       else
01420         uses qh_setdelaunay to project points to paraboloid
01421 */
01422 void qh_projectinput(void) {
01423   int k,i;
01424   int newdim= qh input_dim, newnum= qh num_points;
01425   signed char *project;
01426   int size= (qh input_dim+1)*sizeof(*project);
01427   pointT *newpoints, *coord, *infinity;
01428   realT paraboloid, maxboloid= 0;
01429 
01430   project= (signed char*)qh_memalloc(size);
01431   memset((char*)project, 0, (size_t)size);
01432   for (k=0; k < qh input_dim; k++) {   /* skip Delaunay bound */
01433     if (qh lower_bound[k] == 0 && qh upper_bound[k] == 0) {
01434       project[k]= -1;
01435       newdim--;
01436     }
01437   }
01438   if (qh DELAUNAY) {
01439     project[k]= 1;
01440     newdim++;
01441     if (qh ATinfinity)
01442       newnum++;
01443   }
01444   if (newdim != qh hull_dim) {
01445     qh_fprintf(qh ferr, 6015, "qhull internal error (qh_projectinput): dimension after projection %d != hull_dim %d\n", newdim, qh hull_dim);
01446     qh_errexit(qh_ERRqhull, NULL, NULL);
01447   }
01448   if (!(newpoints=(coordT*)qh_malloc(newnum*newdim*sizeof(coordT)))){
01449     qh_fprintf(qh ferr, 6016, "qhull error: insufficient memory to project %d points\n",
01450            qh num_points);
01451     qh_errexit(qh_ERRmem, NULL, NULL);
01452   }
01453   qh_projectpoints(project, qh input_dim+1, qh first_point,
01454                     qh num_points, qh input_dim, newpoints, newdim);
01455   trace1((qh ferr, 1003, "qh_projectinput: updating lower and upper_bound\n"));
01456   qh_projectpoints(project, qh input_dim+1, qh lower_bound,
01457                     1, qh input_dim+1, qh lower_bound, newdim+1);
01458   qh_projectpoints(project, qh input_dim+1, qh upper_bound,
01459                     1, qh input_dim+1, qh upper_bound, newdim+1);
01460   if (qh HALFspace) {
01461     if (!qh feasible_point) {
01462       qh_fprintf(qh ferr, 6017, "qhull internal error (qh_projectinput): HALFspace defined without qh.feasible_point\n");
01463       qh_errexit(qh_ERRqhull, NULL, NULL);
01464     }
01465     qh_projectpoints(project, qh input_dim, qh feasible_point,
01466                       1, qh input_dim, qh feasible_point, newdim);
01467   }
01468   qh_memfree(project, (qh input_dim+1)*sizeof(*project));
01469   if (qh POINTSmalloc)
01470     qh_free(qh first_point);
01471   qh first_point= newpoints;
01472   qh POINTSmalloc= True;
01473   if (qh DELAUNAY && qh ATinfinity) {
01474     coord= qh first_point;
01475     infinity= qh first_point + qh hull_dim * qh num_points;
01476     for (k=qh hull_dim-1; k--; )
01477       infinity[k]= 0.0;
01478     for (i=qh num_points; i--; ) {
01479       paraboloid= 0.0;
01480       for (k=0; k < qh hull_dim-1; k++) {
01481         paraboloid += *coord * *coord;
01482         infinity[k] += *coord;
01483         coord++;
01484       }
01485       *(coord++)= paraboloid;
01486       maximize_(maxboloid, paraboloid);
01487     }
01488     /* coord == infinity */
01489     for (k=qh hull_dim-1; k--; )
01490       *(coord++) /= qh num_points;
01491     *(coord++)= maxboloid * 1.1;
01492     qh num_points++;
01493     trace0((qh ferr, 9, "qh_projectinput: projected points to paraboloid for Delaunay\n"));
01494   }else if (qh DELAUNAY)  /* !qh ATinfinity */
01495     qh_setdelaunay( qh hull_dim, qh num_points, qh first_point);
01496 } /* projectinput */
01497 
01498 
01499 /*-<a                             href="qh-geom.htm#TOC"
01500   >-------------------------------</a><a name="projectpoints">-</a>
01501 
01502   qh_projectpoints( project, n, points, numpoints, dim, newpoints, newdim )
01503     project points/numpoints/dim to newpoints/newdim
01504     if project[k] == -1
01505       delete dimension k
01506     if project[k] == 1
01507       add dimension k by duplicating previous column
01508     n is size of project
01509 
01510   notes:
01511     newpoints may be points if only adding dimension at end
01512 
01513   design:
01514     check that 'project' and 'newdim' agree
01515     for each dimension
01516       if project == -1
01517         skip dimension
01518       else
01519         determine start of column in newpoints
01520         determine start of column in points
01521           if project == +1, duplicate previous column
01522         copy dimension (column) from points to newpoints
01523 */
01524 void qh_projectpoints(signed char *project, int n, realT *points,
01525         int numpoints, int dim, realT *newpoints, int newdim) {
01526   int testdim= dim, oldk=0, newk=0, i,j=0,k;
01527   realT *newp, *oldp;
01528 
01529   for (k=0; k < n; k++)
01530     testdim += project[k];
01531   if (testdim != newdim) {
01532     qh_fprintf(qh ferr, 6018, "qhull internal error (qh_projectpoints): newdim %d should be %d after projection\n",
01533       newdim, testdim);
01534     qh_errexit(qh_ERRqhull, NULL, NULL);
01535   }
01536   for (j=0; j<n; j++) {
01537     if (project[j] == -1)
01538       oldk++;
01539     else {
01540       newp= newpoints+newk++;
01541       if (project[j] == +1) {
01542         if (oldk >= dim)
01543           continue;
01544         oldp= points+oldk;
01545       }else
01546         oldp= points+oldk++;
01547       for (i=numpoints; i--; ) {
01548         *newp= *oldp;
01549         newp += newdim;
01550         oldp += dim;
01551       }
01552     }
01553     if (oldk >= dim)
01554       break;
01555   }
01556   trace1((qh ferr, 1004, "qh_projectpoints: projected %d points from dim %d to dim %d\n",
01557     numpoints, dim, newdim));
01558 } /* projectpoints */
01559 
01560 
01561 /*-<a                             href="qh-geom.htm#TOC"
01562   >-------------------------------</a><a name="rotateinput">-</a>
01563 
01564   qh_rotateinput( rows )
01565     rotate input using row matrix
01566     input points given by qh first_point, num_points, hull_dim
01567     assumes rows[dim] is a scratch buffer
01568     if qh POINTSmalloc, overwrites input points, else mallocs a new array
01569 
01570   returns:
01571     rotated input
01572     sets qh POINTSmalloc
01573 
01574   design:
01575     see qh_rotatepoints
01576 */
01577 void qh_rotateinput(realT **rows) {
01578 
01579   if (!qh POINTSmalloc) {
01580     qh first_point= qh_copypoints(qh first_point, qh num_points, qh hull_dim);
01581     qh POINTSmalloc= True;
01582   }
01583   qh_rotatepoints(qh first_point, qh num_points, qh hull_dim, rows);
01584 }  /* rotateinput */
01585 
01586 /*-<a                             href="qh-geom.htm#TOC"
01587   >-------------------------------</a><a name="rotatepoints">-</a>
01588 
01589   qh_rotatepoints( points, numpoints, dim, row )
01590     rotate numpoints points by a d-dim row matrix
01591     assumes rows[dim] is a scratch buffer
01592 
01593   returns:
01594     rotated points in place
01595 
01596   design:
01597     for each point
01598       for each coordinate
01599         use row[dim] to compute partial inner product
01600       for each coordinate
01601         rotate by partial inner product
01602 */
01603 void qh_rotatepoints(realT *points, int numpoints, int dim, realT **row) {
01604   realT *point, *rowi, *coord= NULL, sum, *newval;
01605   int i,j,k;
01606 
01607   if (qh IStracing >= 1)
01608     qh_printmatrix(qh ferr, "qh_rotatepoints: rotate points by", row, dim, dim);
01609   for (point= points, j= numpoints; j--; point += dim) {
01610     newval= row[dim];
01611     for (i=0; i < dim; i++) {
01612       rowi= row[i];
01613       coord= point;
01614       for (sum= 0.0, k= dim; k--; )
01615         sum += *rowi++ * *coord++;
01616       *(newval++)= sum;
01617     }
01618     for (k=dim; k--; )
01619       *(--coord)= *(--newval);
01620   }
01621 } /* rotatepoints */
01622 
01623 
01624 /*-<a                             href="qh-geom.htm#TOC"
01625   >-------------------------------</a><a name="scaleinput">-</a>
01626 
01627   qh_scaleinput()
01628     scale input points using qh low_bound/high_bound
01629     input points given by qh first_point, num_points, hull_dim
01630     if qh POINTSmalloc, overwrites input points, else mallocs a new array
01631 
01632   returns:
01633     scales coordinates of points to low_bound[k], high_bound[k]
01634     sets qh POINTSmalloc
01635 
01636   design:
01637     see qh_scalepoints
01638 */
01639 void qh_scaleinput(void) {
01640 
01641   if (!qh POINTSmalloc) {
01642     qh first_point= qh_copypoints(qh first_point, qh num_points, qh hull_dim);
01643     qh POINTSmalloc= True;
01644   }
01645   qh_scalepoints(qh first_point, qh num_points, qh hull_dim,
01646        qh lower_bound, qh upper_bound);
01647 }  /* scaleinput */
01648 
01649 /*-<a                             href="qh-geom.htm#TOC"
01650   >-------------------------------</a><a name="scalelast">-</a>
01651 
01652   qh_scalelast( points, numpoints, dim, low, high, newhigh )
01653     scale last coordinate to [0,m] for Delaunay triangulations
01654     input points given by points, numpoints, dim
01655 
01656   returns:
01657     changes scale of last coordinate from [low, high] to [0, newhigh]
01658     overwrites last coordinate of each point
01659     saves low/high/newhigh in qh.last_low, etc. for qh_setdelaunay()
01660 
01661   notes:
01662     when called by qh_setdelaunay, low/high may not match actual data
01663 
01664   design:
01665     compute scale and shift factors
01666     apply to last coordinate of each point
01667 */
01668 void qh_scalelast(coordT *points, int numpoints, int dim, coordT low,
01669                    coordT high, coordT newhigh) {
01670   realT scale, shift;
01671   coordT *coord;
01672   int i;
01673   boolT nearzero= False;
01674 
01675   trace4((qh ferr, 4013, "qh_scalelast: scale last coordinate from [%2.2g, %2.2g] to [0,%2.2g]\n",
01676     low, high, newhigh));
01677   qh last_low= low;
01678   qh last_high= high;
01679   qh last_newhigh= newhigh;
01680   scale= qh_divzero(newhigh, high - low,
01681                   qh MINdenom_1, &nearzero);
01682   if (nearzero) {
01683     if (qh DELAUNAY)
01684       qh_fprintf(qh ferr, 6019, "qhull input error: can not scale last coordinate.  Input is cocircular\n   or cospherical.   Use option 'Qz' to add a point at infinity.\n");
01685     else
01686       qh_fprintf(qh ferr, 6020, "qhull input error: can not scale last coordinate.  New bounds [0, %2.2g] are too wide for\nexisting bounds [%2.2g, %2.2g] (width %2.2g)\n",
01687                 newhigh, low, high, high-low);
01688     qh_errexit(qh_ERRinput, NULL, NULL);
01689   }
01690   shift= - low * newhigh / (high-low);
01691   coord= points + dim - 1;
01692   for (i=numpoints; i--; coord += dim)
01693     *coord= *coord * scale + shift;
01694 } /* scalelast */
01695 
01696 /*-<a                             href="qh-geom.htm#TOC"
01697   >-------------------------------</a><a name="scalepoints">-</a>
01698 
01699   qh_scalepoints( points, numpoints, dim, newlows, newhighs )
01700     scale points to new lowbound and highbound
01701     retains old bound when newlow= -REALmax or newhigh= +REALmax
01702 
01703   returns:
01704     scaled points
01705     overwrites old points
01706 
01707   design:
01708     for each coordinate
01709       compute current low and high bound
01710       compute scale and shift factors
01711       scale all points
01712       enforce new low and high bound for all points
01713 */
01714 void qh_scalepoints(pointT *points, int numpoints, int dim,
01715         realT *newlows, realT *newhighs) {
01716   int i,k;
01717   realT shift, scale, *coord, low, high, newlow, newhigh, mincoord, maxcoord;
01718   boolT nearzero= False;
01719 
01720   for (k=0; k < dim; k++) {
01721     newhigh= newhighs[k];
01722     newlow= newlows[k];
01723     if (newhigh > REALmax/2 && newlow < -REALmax/2)
01724       continue;
01725     low= REALmax;
01726     high= -REALmax;
01727     for (i=numpoints, coord=points+k; i--; coord += dim) {
01728       minimize_(low, *coord);
01729       maximize_(high, *coord);
01730     }
01731     if (newhigh > REALmax/2)
01732       newhigh= high;
01733     if (newlow < -REALmax/2)
01734       newlow= low;
01735     if (qh DELAUNAY && k == dim-1 && newhigh < newlow) {
01736       qh_fprintf(qh ferr, 6021, "qhull input error: 'Qb%d' or 'QB%d' inverts paraboloid since high bound %.2g < low bound %.2g\n",
01737                k, k, newhigh, newlow);
01738       qh_errexit(qh_ERRinput, NULL, NULL);
01739     }
01740     scale= qh_divzero(newhigh - newlow, high - low,
01741                   qh MINdenom_1, &nearzero);
01742     if (nearzero) {
01743       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",
01744               k, newlow, newhigh, low, high);
01745       qh_errexit(qh_ERRinput, NULL, NULL);
01746     }
01747     shift= (newlow * high - low * newhigh)/(high-low);
01748     coord= points+k;
01749     for (i=numpoints; i--; coord += dim)
01750       *coord= *coord * scale + shift;
01751     coord= points+k;
01752     if (newlow < newhigh) {
01753       mincoord= newlow;
01754       maxcoord= newhigh;
01755     }else {
01756       mincoord= newhigh;
01757       maxcoord= newlow;
01758     }
01759     for (i=numpoints; i--; coord += dim) {
01760       minimize_(*coord, maxcoord);  /* because of roundoff error */
01761       maximize_(*coord, mincoord);
01762     }
01763     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",
01764       k, low, high, newlow, newhigh, numpoints, scale, shift));
01765   }
01766 } /* scalepoints */
01767 
01768 
01769 /*-<a                             href="qh-geom.htm#TOC"
01770   >-------------------------------</a><a name="setdelaunay">-</a>
01771 
01772   qh_setdelaunay( dim, count, points )
01773     project count points to dim-d paraboloid for Delaunay triangulation
01774 
01775     dim is one more than the dimension of the input set
01776     assumes dim is at least 3 (i.e., at least a 2-d Delaunay triangulation)
01777 
01778     points is a dim*count realT array.  The first dim-1 coordinates
01779     are the coordinates of the first input point.  array[dim] is
01780     the first coordinate of the second input point.  array[2*dim] is
01781     the first coordinate of the third input point.
01782 
01783     if qh.last_low defined (i.e., 'Qbb' called qh_scalelast)
01784       calls qh_scalelast to scale the last coordinate the same as the other points
01785 
01786   returns:
01787     for each point
01788       sets point[dim-1] to sum of squares of coordinates
01789     scale points to 'Qbb' if needed
01790 
01791   notes:
01792     to project one point, use
01793       qh_setdelaunay(qh hull_dim, 1, point)
01794 
01795     Do not use options 'Qbk', 'QBk', or 'QbB' since they scale
01796     the coordinates after the original projection.
01797 
01798 */
01799 void qh_setdelaunay(int dim, int count, pointT *points) {
01800   int i, k;
01801   coordT *coordp, coord;
01802   realT paraboloid;
01803 
01804   trace0((qh ferr, 11, "qh_setdelaunay: project %d points to paraboloid for Delaunay triangulation\n", count));
01805   coordp= points;
01806   for (i=0; i < count; i++) {
01807     coord= *coordp++;
01808     paraboloid= coord*coord;
01809     for (k=dim-2; k--; ) {
01810       coord= *coordp++;
01811       paraboloid += coord*coord;
01812     }
01813     *coordp++ = paraboloid;
01814   }
01815   if (qh last_low < REALmax/2)
01816     qh_scalelast(points, count, dim, qh last_low, qh last_high, qh last_newhigh);
01817 } /* setdelaunay */
01818 
01819 
01820 /*-<a                             href="qh-geom.htm#TOC"
01821   >-------------------------------</a><a name="sethalfspace">-</a>
01822 
01823   qh_sethalfspace( dim, coords, nextp, normal, offset, feasible )
01824     set point to dual of halfspace relative to feasible point
01825     halfspace is normal coefficients and offset.
01826 
01827   returns:
01828     false if feasible point is outside of hull (error message already reported)
01829     overwrites coordinates for point at dim coords
01830     nextp= next point (coords)
01831 
01832   design:
01833     compute distance from feasible point to halfspace
01834     divide each normal coefficient by -dist
01835 */
01836 boolT qh_sethalfspace(int dim, coordT *coords, coordT **nextp,
01837          coordT *normal, coordT *offset, coordT *feasible) {
01838   coordT *normp= normal, *feasiblep= feasible, *coordp= coords;
01839   realT dist;
01840   realT r; /*bug fix*/
01841   int k;
01842   boolT zerodiv;
01843 
01844   dist= *offset;
01845   for (k=dim; k--; )
01846     dist += *(normp++) * *(feasiblep++);
01847   if (dist > 0)
01848     goto LABELerroroutside;
01849   normp= normal;
01850   if (dist < -qh MINdenom) {
01851     for (k=dim; k--; )
01852       *(coordp++)= *(normp++) / -dist;
01853   }else {
01854     for (k=dim; k--; ) {
01855       *(coordp++)= qh_divzero(*(normp++), -dist, qh MINdenom_1, &zerodiv);
01856       if (zerodiv)
01857         goto LABELerroroutside;
01858     }
01859   }
01860   *nextp= coordp;
01861   if (qh IStracing >= 4) {
01862     qh_fprintf(qh ferr, 8021, "qh_sethalfspace: halfspace at offset %6.2g to point: ", *offset);
01863     for (k=dim, coordp=coords; k--; ) {
01864       r= *coordp++;
01865       qh_fprintf(qh ferr, 8022, " %6.2g", r);
01866     }
01867     qh_fprintf(qh ferr, 8023, "\n");
01868   }
01869   return True;
01870 LABELerroroutside:
01871   feasiblep= feasible;
01872   normp= normal;
01873   qh_fprintf(qh ferr, 6023, "qhull input error: feasible point is not clearly inside halfspace\nfeasible point: ");
01874   for (k=dim; k--; )
01875     qh_fprintf(qh ferr, 8024, qh_REAL_1, r=*(feasiblep++));
01876   qh_fprintf(qh ferr, 8025, "\n     halfspace: ");
01877   for (k=dim; k--; )
01878     qh_fprintf(qh ferr, 8026, qh_REAL_1, r=*(normp++));
01879   qh_fprintf(qh ferr, 8027, "\n     at offset: ");
01880   qh_fprintf(qh ferr, 8028, qh_REAL_1, *offset);
01881   qh_fprintf(qh ferr, 8029, " and distance: ");
01882   qh_fprintf(qh ferr, 8030, qh_REAL_1, dist);
01883   qh_fprintf(qh ferr, 8031, "\n");
01884   return False;
01885 } /* sethalfspace */
01886 
01887 /*-<a                             href="qh-geom.htm#TOC"
01888   >-------------------------------</a><a name="sethalfspace_all">-</a>
01889 
01890   qh_sethalfspace_all( dim, count, halfspaces, feasible )
01891     generate dual for halfspace intersection with feasible point
01892     array of count halfspaces
01893       each halfspace is normal coefficients followed by offset
01894       the origin is inside the halfspace if the offset is negative
01895 
01896   returns:
01897     malloc'd array of count X dim-1 points
01898 
01899   notes:
01900     call before qh_init_B or qh_initqhull_globals
01901     unused/untested code: please email bradb@shore.net if this works ok for you
01902     If using option 'Fp', also set qh feasible_point. It is a malloc'd array
01903       that is freed by qh_freebuffers.
01904 
01905   design:
01906     see qh_sethalfspace
01907 */
01908 coordT *qh_sethalfspace_all(int dim, int count, coordT *halfspaces, pointT *feasible) {
01909   int i, newdim;
01910   pointT *newpoints;
01911   coordT *coordp, *normalp, *offsetp;
01912 
01913   trace0((qh ferr, 12, "qh_sethalfspace_all: compute dual for halfspace intersection\n"));
01914   newdim= dim - 1;
01915   if (!(newpoints=(coordT*)qh_malloc(count*newdim*sizeof(coordT)))){
01916     qh_fprintf(qh ferr, 6024, "qhull error: insufficient memory to compute dual of %d halfspaces\n",
01917           count);
01918     qh_errexit(qh_ERRmem, NULL, NULL);
01919   }
01920   coordp= newpoints;
01921   normalp= halfspaces;
01922   for (i=0; i < count; i++) {
01923     offsetp= normalp + newdim;
01924     if (!qh_sethalfspace(newdim, coordp, &coordp, normalp, offsetp, feasible)) {
01925       qh_fprintf(qh ferr, 8032, "The halfspace was at index %d\n", i);
01926       qh_errexit(qh_ERRinput, NULL, NULL);
01927     }
01928     normalp= offsetp + 1;
01929   }
01930   return newpoints;
01931 } /* sethalfspace_all */
01932 
01933 
01934 /*-<a                             href="qh-geom.htm#TOC"
01935   >-------------------------------</a><a name="sharpnewfacets">-</a>
01936 
01937   qh_sharpnewfacets()
01938 
01939   returns:
01940     true if could be an acute angle (facets in different quadrants)
01941 
01942   notes:
01943     for qh_findbest
01944 
01945   design:
01946     for all facets on qh.newfacet_list
01947       if two facets are in different quadrants
01948         set issharp
01949 */
01950 boolT qh_sharpnewfacets() {
01951   facetT *facet;
01952   boolT issharp = False;
01953   int *quadrant, k;
01954 
01955   quadrant= (int*)qh_memalloc(qh hull_dim * sizeof(int));
01956   FORALLfacet_(qh newfacet_list) {
01957     if (facet == qh newfacet_list) {
01958       for (k=qh hull_dim; k--; )
01959         quadrant[ k]= (facet->normal[ k] > 0);
01960     }else {
01961       for (k=qh hull_dim; k--; ) {
01962         if (quadrant[ k] != (facet->normal[ k] > 0)) {
01963           issharp= True;
01964           break;
01965         }
01966       }
01967     }
01968     if (issharp)
01969       break;
01970   }
01971   qh_memfree( quadrant, qh hull_dim * sizeof(int));
01972   trace3((qh ferr, 3001, "qh_sharpnewfacets: %d\n", issharp));
01973   return issharp;
01974 } /* sharpnewfacets */
01975 
01976 /*-<a                             href="qh-geom.htm#TOC"
01977   >-------------------------------</a><a name="voronoi_center">-</a>
01978 
01979   qh_voronoi_center( dim, points )
01980     return Voronoi center for a set of points
01981     dim is the orginal dimension of the points
01982     gh.gm_matrix/qh.gm_row are scratch buffers
01983 
01984   returns:
01985     center as a temporary point
01986     if non-simplicial,
01987       returns center for max simplex of points
01988 
01989   notes:
01990     from Bowyer & Woodwark, A Programmer's Geometry, 1983, p. 65
01991 
01992   design:
01993     if non-simplicial
01994       determine max simplex for points
01995     translate point0 of simplex to origin
01996     compute sum of squares of diagonal
01997     compute determinate
01998     compute Voronoi center (see Bowyer & Woodwark)
01999 */
02000 pointT *qh_voronoi_center(int dim, setT *points) {
02001   pointT *point, **pointp, *point0;
02002   pointT *center= (pointT*)qh_memalloc(qh center_size);
02003   setT *simplex;
02004   int i, j, k, size= qh_setsize(points);
02005   coordT *gmcoord;
02006   realT *diffp, sum2, *sum2row, *sum2p, det, factor;
02007   boolT nearzero, infinite;
02008 
02009   if (size == dim+1)
02010     simplex= points;
02011   else if (size < dim+1) {
02012     qh_fprintf(qh ferr, 6025, "qhull internal error (qh_voronoi_center):\n  need at least %d points to construct a Voronoi center\n",
02013              dim+1);
02014     qh_errexit(qh_ERRqhull, NULL, NULL);
02015     simplex= points;  /* never executed -- avoids warning */
02016   }else {
02017     simplex= qh_settemp(dim+1);
02018     qh_maxsimplex(dim, points, NULL, 0, &simplex);
02019   }
02020   point0= SETfirstt_(simplex, pointT);
02021   gmcoord= qh gm_matrix;
02022   for (k=0; k < dim; k++) {
02023     qh gm_row[k]= gmcoord;
02024     FOREACHpoint_(simplex) {
02025       if (point != point0)
02026         *(gmcoord++)= point[k] - point0[k];
02027     }
02028   }
02029   sum2row= gmcoord;
02030   for (i=0; i < dim; i++) {
02031     sum2= 0.0;
02032     for (k=0; k < dim; k++) {
02033       diffp= qh gm_row[k] + i;
02034       sum2 += *diffp * *diffp;
02035     }
02036     *(gmcoord++)= sum2;
02037   }
02038   det= qh_determinant(qh gm_row, dim, &nearzero);
02039   factor= qh_divzero(0.5, det, qh MINdenom, &infinite);
02040   if (infinite) {
02041     for (k=dim; k--; )
02042       center[k]= qh_INFINITE;
02043     if (qh IStracing)
02044       qh_printpoints(qh ferr, "qh_voronoi_center: at infinity for ", simplex);
02045   }else {
02046     for (i=0; i < dim; i++) {
02047       gmcoord= qh gm_matrix;
02048       sum2p= sum2row;
02049       for (k=0; k < dim; k++) {
02050         qh gm_row[k]= gmcoord;
02051         if (k == i) {
02052           for (j=dim; j--; )
02053             *(gmcoord++)= *sum2p++;
02054         }else {
02055           FOREACHpoint_(simplex) {
02056             if (point != point0)
02057               *(gmcoord++)= point[k] - point0[k];
02058           }
02059         }
02060       }
02061       center[i]= qh_determinant(qh gm_row, dim, &nearzero)*factor + point0[i];
02062     }
02063 #ifndef qh_NOtrace
02064     if (qh IStracing >= 3) {
02065       qh_fprintf(qh ferr, 8033, "qh_voronoi_center: det %2.2g factor %2.2g ", det, factor);
02066       qh_printmatrix(qh ferr, "center:", &center, 1, dim);
02067       if (qh IStracing >= 5) {
02068         qh_printpoints(qh ferr, "points", simplex);
02069         FOREACHpoint_(simplex)
02070           qh_fprintf(qh ferr, 8034, "p%d dist %.2g, ", qh_pointid(point),
02071                    qh_pointdist(point, center, dim));
02072         qh_fprintf(qh ferr, 8035, "\n");
02073       }
02074     }
02075 #endif
02076   }
02077   if (simplex != points)
02078     qh_settempfree(&simplex);
02079   return center;
02080 } /* voronoi_center */
02081
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

libqhull
Author(s): Robert Krug
autogenerated on Tue Jun 18 2013 12:38:49