Doxygen Source Code Documentation

Main Page Alphabetical List Data Structures File List Data Fields Globals Search
geom2.c

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

Sections

Personal tools

Navigation

Doxygen Source Code Documentation

geom2.c
