graspit: /opt/ros/diamondback/stacks/graspit_simulator/graspit/graspit

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