00025                                                                      {
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 */

00050                                                                            {
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 */

00078                                                               {
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 */

00125                                                                   {
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 */

00188                            {
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 */

00301                                                                           {
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 */

00345                                                                            {
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 */

00374                                                                   {
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 */

00407                                                                              {
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 */

00454                                    {
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 */

00514                                                                              {
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 */

00587                                         {
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 */

00623                                           {
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 */

00685                                    {
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 */

00755                                     {
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 */

00810                                             {
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 */

00854                                                      {
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 */

00917                            {
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 */

00979                                              {
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 */

01021                                                               {
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 */

01107                          {
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 */

01139                                                                                              {
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 */

01239                                             {
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 */

01259                                                    {
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 */

01286                                        {
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 */

01320                                                                          {
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 */

01361                                                              {
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 */

01386                                                                                    {
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 */

01410                                                            {
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 */

01465                             {
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 */

01568                                                               {
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 */

01620                    {
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 */

01659                              {
01660   realT randr;
01661 
01662   randr= qh_RANDOMint;
01663   return randr * qh RANDOMa + qh RANDOMb;
01664 } /* randomfactor */

01678                                                             {
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 */

01711                                    {
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 */

01737                                                                           {
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 */  

01773                           {
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 */

01803                                                 {
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 */

01849                                          {
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 */    

01933                                                          {
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 */

01971                                                            {
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 */

02042                                                                                        {
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 */

01641                          {
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 */

02092                                                   {
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 */