Skip to content

AFNI/NIfTI Server

Sections
Personal tools
You are here: Home » AFNI » Documentation
 

Doxygen Source Code Documentation

Main Page   Alphabetical List   Data Structures   File List   Data Fields   Globals   Search  

geom.c

Go to the documentation of this file.
00001 /*<html><pre>  -<a                             href="qh-geom.htm"
00002   >-------------------------------</a><a name="TOP">-</a>
00003 
00004    geom.c 
00005    geometric routines of qhull
00006 
00007    see qh-geom.htm and geom.h
00008 
00009    copyright (c) 1993-2001 The Geometry Center        
00010 
00011    infrequent code goes into geom2.c
00012 */
00013    
00014 #include "qhull_a.h"
00015    
00016 /*-<a                             href="qh-geom.htm#TOC"
00017   >-------------------------------</a><a name="backnormal">-</a>
00018   
00019   qh_backnormal( rows, numrow, numcol, sign, normal, nearzero )
00020     given an upper-triangular rows array and a sign,
00021     solve for normal equation x using back substitution over rows U
00022 
00023   returns:
00024      normal= x
00025       
00026      if will not be able to divzero() when normalized (qh.MINdenom_2 and qh.MINdenom_1_2),
00027        if fails on last row
00028          this means that the hyperplane intersects [0,..,1]
00029          sets last coordinate of normal to sign
00030        otherwise
00031          sets tail of normal to [...,sign,0,...], i.e., solves for b= [0...0]
00032          sets nearzero
00033 
00034   notes:
00035      assumes numrow == numcol-1
00036 
00037      see Golub & van Loan 4.4-9 for back substitution
00038 
00039      solves Ux=b where Ax=b and PA=LU
00040      b= [0,...,0,sign or 0]  (sign is either -1 or +1)
00041      last row of A= [0,...,0,1]
00042 
00043      1) Ly=Pb == y=b since P only permutes the 0's of   b
00044      
00045   design:
00046     for each row from end
00047       perform back substitution
00048       if near zero
00049         use qh_divzero for division
00050         if zero divide and not last row
00051           set tail of normal to 0
00052 */
00053 void qh_backnormal (realT **rows, int numrow, int numcol, boolT sign,
00054         coordT *normal, boolT *nearzero) {
00055   int i, j;
00056   coordT *normalp, *normal_tail, *ai, *ak;
00057   realT diagonal;
00058   boolT waszero;
00059   int zerocol= -1;
00060   
00061   normalp= normal + numcol - 1;
00062   *normalp--= (sign ? -1.0 : 1.0);
00063   for(i= numrow; i--; ) {
00064     *normalp= 0.0;
00065     ai= rows[i] + i + 1;
00066     ak= normalp+1;
00067     for(j= i+1; j < numcol; j++)
00068       *normalp -= *ai++ * *ak++;
00069     diagonal= (rows[i])[i];
00070     if (fabs_(diagonal) > qh MINdenom_2)
00071       *(normalp--) /= diagonal;
00072     else {
00073       waszero= False;
00074       *normalp= qh_divzero (*normalp, diagonal, qh MINdenom_1_2, &waszero);
00075       if (waszero) {
00076         zerocol= i;
00077         *(normalp--)= (sign ? -1.0 : 1.0);
00078         for (normal_tail= normalp+2; normal_tail < normal + numcol; normal_tail++)
00079           *normal_tail= 0.0;
00080       }else
00081         normalp--;
00082     }
00083   }
00084   if (zerocol != -1) {
00085     zzinc_(Zback0);
00086     *nearzero= True;
00087     trace4((qh ferr, "qh_backnormal: zero diagonal at column %d.\n", i));
00088     qh_precision ("zero diagonal on back substitution");
00089   }
00090 } /* backnormal */
00091 
00092 /*-<a                             href="qh-geom.htm#TOC"
00093   >-------------------------------</a><a name="distplane">-</a>
00094   
00095   qh_distplane( point, facet, dist )
00096     return distance from point to facet
00097 
00098   returns:
00099     dist
00100     if qh.RANDOMdist, joggles result
00101   
00102   notes:  
00103     dist > 0 if point is above facet (i.e., outside)
00104     does not error (for sortfacets)
00105     
00106   see:
00107     qh_distnorm in geom2.c
00108 */
00109 void qh_distplane (pointT *point, facetT *facet, realT *dist) {
00110   coordT *normal= facet->normal, *coordp, randr;
00111   int k;
00112   
00113   switch(qh hull_dim){
00114   case 2:
00115     *dist= facet->offset + point[0] * normal[0] + point[1] * normal[1];
00116     break;
00117   case 3:
00118     *dist= facet->offset + point[0] * normal[0] + point[1] * normal[1] + point[2] * normal[2];
00119     break;
00120   case 4:
00121     *dist= facet->offset+point[0]*normal[0]+point[1]*normal[1]+point[2]*normal[2]+point[3]*normal[3];
00122     break;
00123   case 5:
00124     *dist= facet->offset+point[0]*normal[0]+point[1]*normal[1]+point[2]*normal[2]+point[3]*normal[3]+point[4]*normal[4];
00125     break;
00126   case 6:
00127     *dist= facet->offset+point[0]*normal[0]+point[1]*normal[1]+point[2]*normal[2]+point[3]*normal[3]+point[4]*normal[4]+point[5]*normal[5];
00128     break;
00129   case 7:  
00130     *dist= facet->offset+point[0]*normal[0]+point[1]*normal[1]+point[2]*normal[2]+point[3]*normal[3]+point[4]*normal[4]+point[5]*normal[5]+point[6]*normal[6];
00131     break;
00132   case 8:
00133     *dist= facet->offset+point[0]*normal[0]+point[1]*normal[1]+point[2]*normal[2]+point[3]*normal[3]+point[4]*normal[4]+point[5]*normal[5]+point[6]*normal[6]+point[7]*normal[7];
00134     break;
00135   default:
00136     *dist= facet->offset;
00137     coordp= point;
00138     for (k= qh hull_dim; k--; )
00139       *dist += *coordp++ * *normal++;
00140     break;
00141   }
00142   zinc_(Zdistplane);
00143   if (!qh RANDOMdist && qh IStracing < 4)
00144     return;
00145   if (qh RANDOMdist) {
00146     randr= qh_RANDOMint;
00147     *dist += (2.0 * randr / qh_RANDOMmax - 1.0) *
00148       qh RANDOMfactor * qh MAXabs_coord;
00149   }
00150   if (qh IStracing >= 4) {
00151     fprintf (qh ferr, "qh_distplane: ");
00152     fprintf (qh ferr, qh_REAL_1, *dist);
00153     fprintf (qh ferr, "from p%d to f%d\n", qh_pointid(point), facet->id);
00154   }
00155   return;
00156 } /* distplane */
00157 
00158 
00159 /*-<a                             href="qh-geom.htm#TOC"
00160   >-------------------------------</a><a name="findbest">-</a>
00161   
00162   qh_findbest( point, startfacet, bestoutside, newfacts, dist, isoutside, numpart )
00163     find facet that is furthest below a point 
00164     searches neighbors of coplanar and most flipped facets
00165     for upperDelaunay facets
00166       searches if coplanar (within searchdist)
00167       returns facet only if !noupper and clearly above
00168 
00169   input:
00170     starts search at 'startfacet' (can not be flipped)
00171     if !bestoutside, stops at qh.MINoutside (DISTround if precise)
00172 
00173   returns:
00174     best facet
00175     dist is distance to facet
00176     isoutside is true if point is outside of facet
00177     numpart counts the number of distance tests
00178 
00179   see also:
00180     qh_findbestnew()
00181     
00182   notes:
00183     uses qh.visit_id, qh.searchset
00184     caller traces the results
00185     after qh_distplane, this and qh_partitionpoint are the most expensive in 3-d
00186       avoid calls to distplane, function calls and real number operations.
00187 
00188   when called by qh_partitionvisible():
00189     indicated by newfacets and isoutside defined
00190     qh.newfacet_list is list of simplicial, new facets
00191     qh_findbestnew set if qh_findbestsharp returns True (to use qh_findbestnew)
00192     qh.bestfacet_notsharp set if qh_findbestsharp returns False
00193     searches horizon of best facet unless "exact" and !bestoutside
00194     searchdist is 2 * DISTround
00195 
00196   when called by qh_check_maxout()
00197     indicated by bestoutside and !newfacets and isoutside == NULL
00198     startfacet must be closest to the point
00199     updates facet->maxoutside for good, visited facets
00200     may return NULL
00201 
00202     searchdist is qh.max_outside + 2 * DISTround
00203       + max( MINvisible('Vn'), MAXcoplanar('Un'));
00204     This setting is a guess.  It must be at least max_outside + 2*DISTround 
00205     because a facet may have a geometric neighbor across a vertex
00206 
00207   when called by findfacet() and check_bestdist()
00208     indicated by !newfacets and isoutside defined
00209     searchdist is same as qh_check_maxout()
00210     returns best facet in neighborhood of given facet
00211       this is best facet overall if dist > -   qh.MAXcoplanar 
00212         or hull has at least a "spherical" curvature
00213 
00214   design:
00215     initialize and test for early exit
00216     repeat while there are facets to search (searchset)
00217       for each neighbor of facet
00218         exit if outside facet found
00219         restart searchset if neighbor is much better
00220         append flipped and nearby facets to searchset
00221     if point is inside and partitioning
00222       test for new facets with a "sharp" intersection
00223       if so, future calls go to qh_findbestsharp
00224     if testhorizon
00225       also test facet's horizon facet
00226 */
00227 facetT *qh_findbest (pointT *point, facetT *startfacet, 
00228                      boolT bestoutside, boolT newfacets, boolT noupper,
00229                      realT *dist, boolT *isoutside, int *numpart) {
00230   realT bestdist= -REALmax/2 /* avoid underflow */, searchdist;
00231   realT cutoff, mincutoff;  /* skip facets that are too far from point */
00232   facetT *facet, *neighbor, **neighborp, *bestfacet= NULL;
00233   int oldtrace= qh IStracing;
00234   int searchsize= 0; /* non-zero if searchset defined */
00235   boolT newbest;
00236   boolT ischeckmax= bestoutside && !newfacets && !isoutside;
00237   boolT ispartition= newfacets && isoutside;
00238   boolT isfindfacet= !newfacets && isoutside;
00239   boolT testhorizon = ispartition && (bestoutside || qh APPROXhull || qh MERGING);
00240 
00241   if (!ischeckmax && !ispartition && !isfindfacet) {
00242     fprintf (qh ferr, "qhull internal error (qh_findbest): unknown combination of arguments\n");
00243     qh_errexit (qh_ERRqhull, startfacet, NULL);
00244   }
00245   if (qh TRACElevel && qh TRACEpoint >= 0 && qh TRACEpoint == qh_pointid (point)) {
00246     qh IStracing= qh TRACElevel;
00247     fprintf (qh ferr, "qh_findbest: point p%d starting at f%d bestoutside? %d newfacets %d\n",
00248              qh TRACEpoint, startfacet->id, bestoutside, newfacets);
00249     fprintf(qh ferr, "  ischeckmax %d ispartition %d isfindfacet %d testhorizon %d\n", 
00250               ischeckmax, ispartition, isfindfacet, testhorizon);
00251     fprintf (qh ferr, "  Last point added to hull was p%d.", qh furthest_id);
00252     fprintf(qh ferr, "  Last merge was #%d.\n", zzval_(Ztotmerge));
00253   }
00254   if (isoutside)
00255     *isoutside= True;
00256 
00257   if (!startfacet->flipped) {
00258     *numpart= 1;           
00259     qh_distplane (point, startfacet, dist);  /* this code is duplicated below */
00260     if (!startfacet->upperdelaunay || (!noupper && *dist >= qh MINoutside)) {
00261       bestdist= *dist;
00262       bestfacet= startfacet;
00263       if (!bestoutside && *dist >= qh MINoutside) 
00264         goto LABELreturn_best;
00265     }
00266 #if qh_MAXoutside
00267     if (ischeckmax && (!qh ONLYgood || startfacet->good) && *dist > startfacet->maxoutside)
00268       startfacet->maxoutside= *dist;
00269 #endif
00270   }
00271   if (ispartition)
00272     searchdist= 2 * qh DISTround;
00273   else 
00274     searchdist= qh max_outside + 2 * qh DISTround
00275                 + fmax_( qh MINvisible, qh MAXcoplanar);
00276   cutoff= bestdist - searchdist;
00277   mincutoff= 0;
00278   if (ischeckmax) {
00279     mincutoff= -(qh DISTround - fmax_(qh MINvisible, qh MAXcoplanar));
00280     minimize_(cutoff, mincutoff);
00281   }
00282   startfacet->visitid= ++qh visit_id;
00283   facet= startfacet;
00284   do {  /* search neighbors of coplanar, upperdelaunay, and flipped facets */
00285     if (True) {
00286  LABELrestart:  /* directed search whenever improvement > searchdist */
00287       newbest= False;
00288       trace4((qh ferr, "qh_findbest: neighbors of f%d, bestdist %2.2g cutoff %2.2g searchdist %2.2g\n", 
00289                   facet->id, bestdist, cutoff, searchdist));
00290       FOREACHneighbor_(facet) {
00291         if (ispartition && !neighbor->newfacet)
00292           continue;
00293         if (!neighbor->flipped) {
00294           if (neighbor->visitid == qh visit_id)
00295             continue;
00296           neighbor->visitid= qh visit_id;
00297           (*numpart)++;
00298           qh_distplane (point, neighbor, dist);
00299           if (!bestoutside && *dist >= qh MINoutside 
00300           && (!noupper || !facet->upperdelaunay)) {
00301             bestfacet= neighbor;
00302             goto LABELreturn_best;
00303           }
00304 #if qh_MAXoutside
00305           if (ischeckmax) {
00306             if ((!qh ONLYgood || neighbor->good) 
00307             && *dist > neighbor->maxoutside)
00308               neighbor->maxoutside= *dist;
00309             else if (bestfacet && *dist < cutoff)
00310               continue;
00311           }else 
00312 #endif      /* dangling else! */
00313           if (bestfacet && *dist < cutoff)
00314             continue;
00315           if (*dist > bestdist) {
00316             if (!neighbor->upperdelaunay 
00317             || (bestoutside && !noupper && *dist >= qh MINoutside)) {
00318               if (ischeckmax && qh_MAXoutside) {
00319                 bestdist= *dist;
00320                 bestfacet= neighbor;
00321                 cutoff= bestdist - searchdist;
00322                 minimize_(cutoff, mincutoff);
00323               }else if (*dist > bestdist + searchdist) {
00324                 bestdist= *dist;
00325                 bestfacet= neighbor;
00326                 cutoff= bestdist - searchdist;
00327                 searchsize= 0;
00328                 facet= neighbor;
00329                 if (newbest) /* newbest may be coplanar with facet */
00330                   facet->visitid= ++qh visit_id;
00331                 goto LABELrestart; /* repeat with a new facet */
00332               }else {
00333                 bestdist= *dist;
00334                 bestfacet= neighbor;
00335                 cutoff= bestdist - searchdist;
00336               }
00337               newbest= True;
00338             }
00339           }
00340         }
00341         if (!searchsize++) {
00342           SETfirst_(qh searchset) = neighbor;
00343           qh_settruncate (qh searchset, 1);
00344         }else
00345           qh_setappend (&qh searchset, neighbor);
00346       } /* FOREACHneighbor */
00347     } /* while (restart) */
00348   }while
00349     (searchsize && (facet= (facetT*)qh_setdellast (qh searchset)));
00350   if (!ischeckmax) {
00351     if (!bestfacet) {
00352       fprintf (qh ferr, "qh_findbest: point p%d starting at f%d bestoutside? %d newfacets %d\n",
00353                qh TRACEpoint, startfacet->id, bestoutside, newfacets);
00354       fprintf(qh ferr, "\n\
00355 qh_findbest: all neighbors of facet %d are flipped or upper Delaunay.\n\
00356 Please report this error to qhull_bug@geom.umn.edu with the input and all of the output.\n",
00357          startfacet->id);
00358       qh FORCEoutput= True;
00359       qh_errexit (qh_ERRqhull, startfacet, NULL);
00360     }
00361     if (ispartition && !qh findbest_notsharp && bestdist < - qh DISTround) {
00362       if (qh_findbestsharp ( point, &bestfacet, &bestdist, numpart)) 
00363         qh findbestnew= True;
00364       else
00365         qh findbest_notsharp= True;
00366     }
00367     if (testhorizon) {
00368       facet= SETfirstt_(bestfacet->neighbors, facetT);
00369       trace4((qh ferr, "qh_findbest: horizon facet f%d\n", facet->id));
00370       (*numpart)++;
00371       qh_distplane (point, facet, dist);
00372       if (*dist > bestdist 
00373       && (!facet->upperdelaunay || (!noupper && *dist >= qh MINoutside))) {
00374         bestdist= *dist;
00375         bestfacet= facet;
00376       }
00377     }
00378   }
00379   *dist= bestdist;
00380   if (isoutside && bestdist < qh MINoutside)
00381     *isoutside= False;
00382 LABELreturn_best:
00383   qh IStracing= oldtrace;
00384   return bestfacet;
00385 }  /* findbest */
00386 
00387 
00388 /*-<a                             href="qh-geom.htm#TOC"
00389   >-------------------------------</a><a name="findbestnew">-</a>
00390   
00391   qh_findbestnew( point, startfacet, dist, isoutside, numpart )
00392     find best newfacet for point
00393     searches new facets starting at startfacet
00394 
00395   returns:
00396     dist is distance to facet
00397     isoutside is true if point is outside of facet
00398     numpart is number of distance tests
00399 
00400   notes:
00401     if qh.BESToutside or !isoutside
00402        stops at furthest facet
00403     if qh.MERGING 
00404        stops when distance > qh_DISToutside (max(4*MINoutside, 2*max_outside))
00405     else 
00406        stops when distance > MINoutside (DISTround in precise case)
00407     searches newfacets then searchs neighbors of best facet.
00408     avoids upperdelaunay facet unless none other or (isoutside and outside)
00409 
00410     uses visit_id and seen flags
00411     caller traces the results
00412   
00413   see also:
00414     qh_partitionall() and qh_findbest()
00415 
00416   design:
00417     for each new facet
00418       test distance from point to facet
00419       select best facet
00420     test each horizon facet of best facet
00421     return best facet 
00422 */
00423 facetT *qh_findbestnew (pointT *point, facetT *startfacet,
00424            realT *dist, boolT *isoutside, int *numpart) {
00425   realT bestdist= -REALmax, bestdist2= -REALmax;
00426   facetT *neighbor, **neighborp, *bestfacet= NULL, *newfacet, *facet;
00427   facetT *bestfacet2= NULL;
00428   int oldtrace= qh IStracing, i;
00429   realT distoutside;
00430 
00431   if (!startfacet) {
00432     if (qh MERGING)
00433       fprintf(qh ferr, "qhull precision error (qh_findbestnew): merging has formed and deleted an independent cycle of facets.  Can not continue.\n");
00434     else
00435       fprintf(qh ferr, "qhull internal error (qh_findbestnew): no new facets for point p%d\n",
00436               qh furthest_id);      
00437     qh_errexit (qh_ERRqhull, NULL, NULL);
00438   }
00439   if (qh BESToutside || !isoutside)
00440     distoutside= REALmax;
00441   else if (qh MERGING)
00442     distoutside= qh_DISToutside; /* defined in user.h */
00443   else
00444     distoutside= qh MINoutside;
00445   if (qh TRACElevel && qh TRACEpoint >= 0 && qh TRACEpoint == qh_pointid (point)) {
00446     qh IStracing= qh TRACElevel;
00447     fprintf(qh ferr, "qh_findbestnew: point p%d facet f%d. Stop if dist > %2.2g\n",
00448              qh TRACEpoint, startfacet->id, distoutside);
00449     fprintf(qh ferr, "  Last point added to hull was p%d.", qh furthest_id);
00450     fprintf(qh ferr, "  Last merge was #%d.\n", zzval_(Ztotmerge));
00451   }
00452   if (isoutside)
00453     *isoutside= True;
00454   *numpart= 0;
00455 
00456   /* visit all new facets starting with startfacet */
00457   for (i= 0, facet= startfacet; i < 2; i++, facet= qh newfacet_list) {
00458     FORALLfacet_(facet) {
00459       if (facet == startfacet && i)
00460         break;
00461       qh_distplane (point, facet, dist);
00462       (*numpart)++;
00463       if (facet->upperdelaunay) {
00464         if (*dist > bestdist2) {
00465           bestdist2= *dist;
00466           bestfacet2= facet;
00467           if (*dist >= distoutside) {
00468             bestfacet= facet;
00469             goto LABELreturn_bestnew;
00470           }
00471         }
00472       }else if (*dist > bestdist) {
00473         bestdist= *dist;
00474         bestfacet= facet;
00475         if (*dist >= distoutside) 
00476           goto LABELreturn_bestnew;
00477       }
00478     }
00479   }
00480   newfacet= bestfacet ? bestfacet : bestfacet2;
00481         /* !bestfacet only occurs if 'd' creates incorrect upper-delaunay facets */
00482   FOREACHneighbor_(newfacet) {
00483     if (!neighbor->newfacet) {
00484       qh_distplane (point, neighbor, dist);
00485       (*numpart)++;
00486       if (neighbor->upperdelaunay) {
00487         if (*dist > bestdist2) {
00488           bestdist2= *dist;
00489           bestfacet2= neighbor;
00490         }
00491       }else if (*dist > bestdist) {
00492         bestdist= *dist;
00493         bestfacet= neighbor;
00494       }
00495     }
00496   }
00497   if (!bestfacet  
00498   || (isoutside && bestdist2 >= qh MINoutside && bestdist2 > bestdist)) {
00499     *dist= bestdist2;
00500     bestfacet= bestfacet2;
00501   }else 
00502     *dist= bestdist;
00503   if (isoutside && *dist < qh MINoutside)
00504     *isoutside= False;
00505 LABELreturn_bestnew:
00506   qh IStracing= oldtrace;
00507   return bestfacet;
00508 }  /* findbestnew */
00509 
00510 
00511 /*-<a                             href="qh-geom.htm#TOC"
00512   >-------------------------------</a><a name="gausselim">-</a>
00513   
00514   qh_gausselim( rows, numrow, numcol, sign )
00515     Gaussian elimination with partial pivoting
00516 
00517   returns:
00518     rows is upper triangular (includes row exchanges)
00519     flips sign for each row exchange
00520     sets nearzero if pivot[k] < qh.NEARzero[k], else clears it
00521 
00522   notes:
00523     if nearzero, the determinant's sign may be incorrect.
00524     assumes numrow <= numcol
00525 
00526   design:
00527     for each row
00528       determine pivot and exchange rows if necessary
00529       test for near zero
00530       perform gaussian elimination step
00531 */
00532 void qh_gausselim(realT **rows, int numrow, int numcol, boolT *sign, boolT *nearzero) {
00533   realT *ai, *ak, *rowp, *pivotrow;
00534   realT n, pivot, pivot_abs= 0.0, temp;
00535   int i, j, k, pivoti, flip=0;
00536   
00537   *nearzero= False;
00538   for(k= 0; k < numrow; k++) {
00539     pivot_abs= fabs_((rows[k])[k]);
00540     pivoti= k;
00541     for(i= k+1; i < numrow; i++) {
00542       if ((temp= fabs_((rows[i])[k])) > pivot_abs) {
00543         pivot_abs= temp;
00544         pivoti= i;
00545       }
00546     }
00547     if (pivoti != k) {
00548       rowp= rows[pivoti]; 
00549       rows[pivoti]= rows[k]; 
00550       rows[k]= rowp; 
00551       *sign ^= 1;
00552       flip ^= 1;
00553     }
00554     if (pivot_abs <= qh NEARzero[k]) {
00555       *nearzero= True;
00556       if (pivot_abs == 0.0) {   /* remainder of column == 0 */
00557         if (qh IStracing >= 4) {
00558           fprintf (qh ferr, "qh_gausselim: 0 pivot at column %d. (%2.2g < %2.2g)\n", k, pivot_abs, qh DISTround);
00559           qh_printmatrix (qh ferr, "Matrix:", rows, numrow, numcol);
00560         }
00561         zzinc_(Zgauss0);
00562         qh_precision ("zero pivot for Gaussian elimination");
00563         goto LABELnextcol;
00564       }
00565     }
00566     pivotrow= rows[k] + k;
00567     pivot= *pivotrow++;  /* signed value of pivot, and remainder of row */
00568     for(i= k+1; i < numrow; i++) {
00569       ai= rows[i] + k;
00570       ak= pivotrow;
00571       n= (*ai++)/pivot;   /* divzero() not needed since |pivot| >= |*ai| */
00572       for(j= numcol - (k+1); j--; )
00573         *ai++ -= n * *ak++;
00574     }
00575   LABELnextcol:
00576     ;
00577   }
00578   wmin_(Wmindenom, pivot_abs);  /* last pivot element */
00579   if (qh IStracing >= 5)
00580     qh_printmatrix (qh ferr, "qh_gausselem: result", rows, numrow, numcol);
00581 } /* gausselim */
00582 
00583 
00584 /*-<a                             href="qh-geom.htm#TOC"
00585   >-------------------------------</a><a name="getangle">-</a>
00586   
00587   qh_getangle( vect1, vect2 )
00588     returns the dot product of two, DIM3 vectors
00589     if qh.RANDOMdist, joggles result
00590 
00591   notes:
00592     the angle may be > 1.0 or < -1.0 because of roundoff errors
00593 
00594 */
00595 realT qh_getangle(pointT *vect1, pointT *vect2) {
00596   realT angle= 0, randr;
00597   int k;
00598 
00599   for(k= qh hull_dim; k--; )
00600     angle += *vect1++ * *vect2++;
00601   if (qh RANDOMdist) {
00602     randr= qh_RANDOMint;
00603     angle += (2.0 * randr / qh_RANDOMmax - 1.0) *
00604       qh RANDOMfactor;
00605   }
00606   trace4((qh ferr, "qh_getangle: %2.2g\n", angle));
00607   return(angle);
00608 } /* getangle */
00609 
00610 
00611 /*-<a                             href="qh-geom.htm#TOC"
00612   >-------------------------------</a><a name="getcenter">-</a>
00613   
00614   qh_getcenter( vertices )
00615     returns arithmetic center of a set of vertices as a new point
00616 
00617   notes:
00618     allocates point array for center
00619 */
00620 pointT *qh_getcenter(setT *vertices) {
00621   int k;
00622   pointT *center, *coord;
00623   vertexT *vertex, **vertexp;
00624   int count= qh_setsize(vertices);
00625 
00626   if (count < 2) {
00627     fprintf (qh ferr, "qhull internal error (qh_getcenter): not defined for %d points\n", count);
00628     qh_errexit (qh_ERRqhull, NULL, NULL);
00629   }
00630   center= (pointT *)qh_memalloc(qh normal_size);
00631   for (k=0; k < qh hull_dim; k++) {
00632     coord= center+k;
00633     *coord= 0.0;
00634     FOREACHvertex_(vertices)
00635       *coord += vertex->point[k];
00636     *coord /= count;
00637   }
00638   return(center);
00639 } /* getcenter */
00640 
00641 
00642 /*-<a                             href="qh-geom.htm#TOC"
00643   >-------------------------------</a><a name="getcentrum">-</a>
00644   
00645   qh_getcentrum( facet )
00646     returns the centrum for a facet as a new point
00647 
00648   notes:
00649     allocates the centrum
00650 */
00651 pointT *qh_getcentrum(facetT *facet) {
00652   realT dist;
00653   pointT *centrum, *point;
00654 
00655   point= qh_getcenter(facet->vertices);
00656   zzinc_(Zcentrumtests);
00657   qh_distplane (point, facet, &dist);
00658   centrum= qh_projectpoint(point, facet, dist);
00659   qh_memfree(point, qh normal_size);
00660   trace4((qh ferr, "qh_getcentrum: for f%d, %d vertices dist= %2.2g\n",
00661           facet->id, qh_setsize(facet->vertices), dist));
00662   return centrum;
00663 } /* getcentrum */
00664 
00665 
00666 /*-<a                             href="qh-geom.htm#TOC"
00667   >-------------------------------</a><a name="getdistance">-</a>
00668   
00669   qh_getdistance( facet, neighbor, mindist, maxdist )
00670     returns the maxdist and mindist distance of any vertex from neighbor
00671 
00672   returns:
00673     the max absolute value
00674 
00675   design:
00676     for each vertex of facet that is not in neighbor
00677       test the distance from vertex to neighbor
00678 */
00679 realT qh_getdistance(facetT *facet, facetT *neighbor, realT *mindist, realT *maxdist) {
00680   vertexT *vertex, **vertexp;
00681   realT dist, maxd, mind;
00682   
00683   FOREACHvertex_(facet->vertices)
00684     vertex->seen= False;
00685   FOREACHvertex_(neighbor->vertices)
00686     vertex->seen= True;
00687   mind= 0.0;
00688   maxd= 0.0;
00689   FOREACHvertex_(facet->vertices) {
00690     if (!vertex->seen) {
00691       zzinc_(Zbestdist);
00692       qh_distplane(vertex->point, neighbor, &dist);
00693       if (dist < mind)
00694         mind= dist;
00695       else if (dist > maxd)
00696         maxd= dist;
00697     }
00698   }
00699   *mindist= mind;
00700   *maxdist= maxd;
00701   mind= -mind;
00702   if (maxd > mind)
00703     return maxd;
00704   else
00705     return mind;
00706 } /* getdistance */
00707 
00708 
00709 /*-<a                             href="qh-geom.htm#TOC"
00710   >-------------------------------</a><a name="normalize">-</a>
00711 
00712   qh_normalize( normal, dim, toporient )
00713     normalize a vector and report if too small
00714     does not use min norm
00715   
00716   see:
00717     qh_normalize2
00718 */
00719 void qh_normalize (coordT *normal, int dim, boolT toporient) {
00720   qh_normalize2( normal, dim, toporient, NULL, NULL);
00721 } /* normalize */
00722 
00723 /*-<a                             href="qh-geom.htm#TOC"
00724   >-------------------------------</a><a name="normalize2">-</a>
00725   
00726   qh_normalize2( normal, dim, toporient, minnorm, ismin )
00727     normalize a vector and report if too small
00728     qh.MINdenom/MINdenom1 are the upper limits for divide overflow
00729 
00730   returns:
00731     normalized vector
00732     flips sign if !toporient
00733     if minnorm non-NULL, 
00734       sets ismin if normal < minnorm
00735 
00736   notes:
00737     if zero norm
00738        sets all elements to sqrt(1.0/dim)
00739     if divide by zero (divzero ())
00740        sets largest element to   +/-1
00741        bumps Znearlysingular
00742       
00743   design:
00744     computes norm
00745     test for minnorm
00746     if not near zero
00747       normalizes normal
00748     else if zero norm
00749       sets normal to standard value
00750     else
00751       uses qh_divzero to normalize
00752       if nearzero
00753         sets norm to direction of maximum value
00754 */
00755 void qh_normalize2 (coordT *normal, int dim, boolT toporient, 
00756             realT *minnorm, boolT *ismin) {
00757   int k;
00758   realT *colp, *maxp, norm= 0, temp, *norm1, *norm2, *norm3;
00759   boolT zerodiv;
00760 
00761   norm1= normal+1;
00762   norm2= normal+2;
00763   norm3= normal+3;
00764   if (dim == 2)
00765     norm= sqrt((*normal)*(*normal) + (*norm1)*(*norm1));
00766   else if (dim == 3)
00767     norm= sqrt((*normal)*(*normal) + (*norm1)*(*norm1) + (*norm2)*(*norm2));
00768   else if (dim == 4) {
00769     norm= sqrt((*normal)*(*normal) + (*norm1)*(*norm1) + (*norm2)*(*norm2) 
00770                + (*norm3)*(*norm3));
00771   }else if (dim > 4) {
00772     norm= (*normal)*(*normal) + (*norm1)*(*norm1) + (*norm2)*(*norm2) 
00773                + (*norm3)*(*norm3);
00774     for (k= dim-4, colp= normal+4; k--; colp++)
00775       norm += (*colp) * (*colp);
00776     norm= sqrt(norm);
00777   }
00778   if (minnorm) {
00779     if (norm < *minnorm) 
00780       *ismin= True;
00781     else
00782       *ismin= False;
00783   }
00784   wmin_(Wmindenom, norm);
00785   if (norm > qh MINdenom) {
00786     if (!toporient)
00787       norm= -norm;
00788     *normal /= norm;
00789     *norm1 /= norm;
00790     if (dim == 2)
00791       ; /* all done */
00792     else if (dim == 3)
00793       *norm2 /= norm;
00794     else if (dim == 4) {
00795       *norm2 /= norm;
00796       *norm3 /= norm;
00797     }else if (dim >4) {
00798       *norm2 /= norm;
00799       *norm3 /= norm;
00800       for (k= dim-4, colp= normal+4; k--; )
00801         *colp++ /= norm;
00802     }
00803   }else if (norm == 0.0) {
00804     temp= sqrt (1.0/dim);
00805     for (k= dim, colp= normal; k--; )
00806       *colp++ = temp;
00807   }else {
00808     if (!toporient)
00809       norm= -norm;
00810     for (k= dim, colp= normal; k--; colp++) { /* k used below */
00811       temp= qh_divzero (*colp, norm, qh MINdenom_1, &zerodiv);
00812       if (!zerodiv)
00813         *colp= temp;
00814       else {
00815         maxp= qh_maxabsval(normal, dim);
00816         temp= ((*maxp * norm >= 0.0) ? 1.0 : -1.0);
00817         for (k= dim, colp= normal; k--; colp++)
00818           *colp= 0.0;
00819         *maxp= temp;
00820         zzinc_(Znearlysingular);
00821         trace0((qh ferr, "qh_normalize: norm=%2.2g too small during p%d\n", 
00822                norm, qh furthest_id));
00823         return;
00824       }
00825     }
00826   }
00827 } /* normalize */
00828 
00829 
00830 /*-<a                             href="qh-geom.htm#TOC"
00831   >-------------------------------</a><a name="projectpoint">-</a>
00832   
00833   qh_projectpoint( point, facet, dist )
00834     project point onto a facet by dist
00835 
00836   returns:
00837     returns a new point
00838     
00839   notes:
00840     if dist= distplane(point,facet)
00841       this projects point to hyperplane
00842     assumes qh_memfree_() is valid for normal_size
00843 */
00844 pointT *qh_projectpoint(pointT *point, facetT *facet, realT dist) {
00845   pointT *newpoint, *np, *normal;
00846   int normsize= qh normal_size,k;
00847   void **freelistp; /* used !qh_NOmem */
00848   
00849   qh_memalloc_(normsize, freelistp, newpoint, pointT);
00850   np= newpoint;
00851   normal= facet->normal;
00852   for(k= qh hull_dim; k--; )
00853     *(np++)= *point++ - dist * *normal++;
00854   return(newpoint);
00855 } /* projectpoint */
00856 
00857   
00858 /*-<a                             href="qh-geom.htm#TOC"
00859   >-------------------------------</a><a name="setfacetplane">-</a>
00860   
00861   qh_setfacetplane( facet )
00862     sets the hyperplane for a facet
00863     if qh.RANDOMdist, joggles hyperplane
00864 
00865   notes:
00866     uses global buffers qh.gm_matrix and qh.gm_row
00867     overwrites facet->normal if already defined
00868     updates Wnewvertex if PRINTstatistics
00869     sets facet->upperdelaunay if upper envelope of Delaunay triangulation
00870 
00871   design:
00872     copy vertex coordinates to qh.gm_matrix/gm_row
00873     compute determinate
00874     if nearzero
00875       recompute determinate with gaussian elimination
00876       if nearzero
00877         force outside orientation by testing interior point
00878 */
00879 void qh_setfacetplane(facetT *facet) {
00880   pointT *point;
00881   vertexT *vertex, **vertexp;
00882   int k,i, normsize= qh normal_size, oldtrace= 0;
00883   realT dist;
00884   void **freelistp; /* used !qh_NOmem */
00885   coordT *coord, *gmcoord;
00886   pointT *point0= SETfirstt_(facet->vertices, vertexT)->point;
00887   boolT nearzero= False;
00888 
00889   zzinc_(Zsetplane);
00890   if (!facet->normal)
00891     qh_memalloc_(normsize, freelistp, facet->normal, coordT);
00892   if (facet == qh tracefacet) {
00893     oldtrace= qh IStracing;
00894     qh IStracing= 5;
00895     fprintf (qh ferr, "qh_setfacetplane: facet f%d created.\n", facet->id);
00896     fprintf (qh ferr, "  Last point added to hull was p%d.", qh furthest_id);
00897     if (zzval_(Ztotmerge))
00898       fprintf(qh ferr, "  Last merge was #%d.", zzval_(Ztotmerge));
00899     fprintf (qh ferr, "\n\nCurrent summary is:\n");
00900       qh_printsummary (qh ferr);
00901   }
00902   if (qh hull_dim <= 4) {
00903     i= 0;
00904     if (qh RANDOMdist) {
00905       gmcoord= qh gm_matrix;
00906       FOREACHvertex_(facet->vertices) {
00907         qh gm_row[i++]= gmcoord;
00908         coord= vertex->point;
00909         for (k= qh hull_dim; k--; )
00910           *(gmcoord++)= *coord++ * qh_randomfactor();
00911       }   
00912     }else {
00913       FOREACHvertex_(facet->vertices)
00914        qh gm_row[i++]= vertex->point;
00915     }
00916     qh_sethyperplane_det(qh hull_dim, qh gm_row, point0, facet->toporient,
00917                 facet->normal, &facet->offset, &nearzero);
00918   }
00919   if (qh hull_dim > 4 || nearzero) {
00920     i= 0;
00921     gmcoord= qh gm_matrix;
00922     FOREACHvertex_(facet->vertices) {
00923       if (vertex->point != point0) {
00924         qh gm_row[i++]= gmcoord;
00925         coord= vertex->point;
00926         point= point0;
00927         for(k= qh hull_dim; k--; )
00928           *(gmcoord++)= *coord++ - *point++;
00929       }
00930     }
00931     qh gm_row[i]= gmcoord;  /* for areasimplex */
00932     if (qh RANDOMdist) {
00933       gmcoord= qh gm_matrix;
00934       for (i= qh hull_dim-1; i--; ) {
00935         for (k= qh hull_dim; k--; )
00936           *(gmcoord++) *= qh_randomfactor();
00937       }
00938     }
00939     qh_sethyperplane_gauss(qh hull_dim, qh gm_row, point0, facet->toporient,
00940                 facet->normal, &facet->offset, &nearzero);
00941     if (nearzero) { 
00942       if (qh_orientoutside (facet)) {
00943         trace0((qh ferr, "qh_setfacetplane: flipped orientation after testing interior_point during p%d\n", qh furthest_id));
00944       /* this is part of using Gaussian Elimination.  For example in 5-d
00945            1 1 1 1 0
00946            1 1 1 1 1
00947            0 0 0 1 0
00948            0 1 0 0 0
00949            1 0 0 0 0
00950            norm= 0.38 0.38 -0.76 0.38 0
00951          has a determinate of 1, but g.e. after subtracting pt. 0 has
00952          0's in the diagonal, even with full pivoting.  It does work
00953          if you subtract pt. 4 instead. */
00954       }
00955     }
00956   }
00957   facet->upperdelaunay= False;
00958   if (qh DELAUNAY) {
00959     if (qh UPPERdelaunay) {     /* matches qh.lower_threshold in qh_initbuild */
00960       if (facet->normal[qh hull_dim -1] >= qh ANGLEround * qh_ZEROdelaunay)
00961         facet->upperdelaunay= True;
00962     }else {
00963       if (facet->normal[qh hull_dim -1] > -qh ANGLEround * qh_ZEROdelaunay)
00964         facet->upperdelaunay= True;
00965     }
00966   }
00967   if (qh PRINTstatistics || qh IStracing || qh TRACElevel || qh JOGGLEmax < REALmax) {
00968     qh old_randomdist= qh RANDOMdist;
00969     qh RANDOMdist= False;
00970     FOREACHvertex_(facet->vertices) {
00971       if (vertex->point != point0) {
00972         boolT istrace= False;
00973         zinc_(Zdiststat);
00974         qh_distplane(vertex->point, facet, &dist);
00975         dist= fabs_(dist);
00976         zinc_(Znewvertex);
00977         wadd_(Wnewvertex, dist);
00978         if (dist > wwval_(Wnewvertexmax)) {
00979           wwval_(Wnewvertexmax)= dist;
00980           if (dist > qh max_outside) {
00981             qh max_outside= dist;  /* used by qh_maxouter() */
00982             if (dist > qh TRACEdist) 
00983               istrace= True;
00984           }
00985         }else if (-dist > qh TRACEdist)
00986           istrace= True;
00987         if (istrace) {
00988           fprintf (qh ferr, "qh_setfacetplane: ====== vertex p%d (v%d) increases max_outside to %2.2g for new facet f%d last p%d\n",
00989                 qh_pointid(vertex->point), vertex->id, dist, facet->id, qh furthest_id);
00990           qh_errprint ("DISTANT", facet, NULL, NULL, NULL);
00991         }
00992       }
00993     }
00994     qh RANDOMdist= qh old_randomdist;
00995   }
00996   if (qh IStracing >= 3) {
00997     fprintf (qh ferr, "qh_setfacetplane: f%d offset %2.2g normal: ",
00998              facet->id, facet->offset);
00999     for (k=0; k < qh hull_dim; k++)
01000       fprintf (qh ferr, "%2.2g ", facet->normal[k]);
01001     fprintf (qh ferr, "\n");
01002   }
01003   if (facet == qh tracefacet)
01004     qh IStracing= oldtrace;
01005 } /* setfacetplane */
01006 
01007 
01008 /*-<a                             href="qh-geom.htm#TOC"
01009   >-------------------------------</a><a name="sethyperplane_det">-</a>
01010   
01011   qh_sethyperplane_det( dim, rows, point0, toporient, normal, offset, nearzero )
01012     given dim X dim array indexed by rows[], one row per point, 
01013         toporient (flips all signs),
01014         and point0 (any row)
01015     set normalized hyperplane equation from oriented simplex
01016 
01017   returns:
01018     normal (normalized)
01019     offset (places point0 on the hyperplane)
01020     sets nearzero if hyperplane not through points
01021 
01022   notes:
01023     only defined for dim == 2..4
01024     rows[] is not modified
01025     solves det(P-V_0, V_n-V_0, ..., V_1-V_0)=0, i.e. every point is on hyperplane
01026     see Bower & Woodworth, A programmer's geometry, Butterworths 1983.
01027 
01028   derivation of 3-d minnorm
01029     Goal: all vertices V_i within qh.one_merge of hyperplane
01030     Plan: exactly translate the facet so that V_0 is the origin
01031           exactly rotate the facet so that V_1 is on the x-axis and y_2=0.
01032           exactly rotate the effective perturbation to only effect n_0
01033              this introduces a factor of sqrt(3)
01034     n_0 = ((y_2-y_0)*(z_1-z_0) - (z_2-z_0)*(y_1-y_0)) / norm
01035     Let M_d be the max coordinate difference
01036     Let M_a be the greater of M_d and the max abs. coordinate
01037     Let u be machine roundoff and distround be max error for distance computation
01038     The max error for n_0 is sqrt(3) u M_a M_d / norm.  n_1 is approx. 1 and n_2 is approx. 0
01039     The max error for distance of V_1 is sqrt(3) u M_a M_d M_d / norm.  Offset=0 at origin
01040     Then minnorm = 1.8 u M_a M_d M_d / qh.ONEmerge
01041     Note that qh.one_merge is approx. 45.5 u M_a and norm is usually about M_d M_d
01042 
01043   derivation of 4-d minnorm
01044     same as above except rotate the facet so that V_1 on x-axis and w_2, y_3, w_3=0
01045      [if two vertices fixed on x-axis, can rotate the other two in yzw.]
01046     n_0 = det3_(...) = y_2 det2_(z_1, w_1, z_3, w_3) = - y_2 w_1 z_3
01047      [all other terms contain at least two factors nearly zero.]
01048     The max error for n_0 is sqrt(4) u M_a M_d M_d / norm
01049     Then minnorm = 2 u M_a M_d M_d M_d / qh.ONEmerge
01050     Note that qh.one_merge is approx. 82 u M_a and norm is usually about M_d M_d M_d
01051 */
01052 void qh_sethyperplane_det (int dim, coordT **rows, coordT *point0, 
01053           boolT toporient, coordT *normal, realT *offset, boolT *nearzero) {
01054   realT maxround, dist;
01055   int i;
01056   pointT *point;
01057 
01058 
01059   if (dim == 2) {
01060     normal[0]= dY(1,0);
01061     normal[1]= dX(0,1);
01062     qh_normalize2 (normal, dim, toporient, NULL, NULL);
01063     *offset= -(point0[0]*normal[0]+point0[1]*normal[1]);
01064     *nearzero= False;  /* since nearzero norm => incident points */
01065   }else if (dim == 3) {
01066     normal[0]= det2_(dY(2,0), dZ(2,0),
01067                      dY(1,0), dZ(1,0));
01068     normal[1]= det2_(dX(1,0), dZ(1,0),
01069                      dX(2,0), dZ(2,0));
01070     normal[2]= det2_(dX(2,0), dY(2,0),
01071                      dX(1,0), dY(1,0));
01072     qh_normalize2 (normal, dim, toporient, NULL, NULL);
01073     *offset= -(point0[0]*normal[0] + point0[1]*normal[1]
01074                + point0[2]*normal[2]);
01075     maxround= qh DISTround;
01076     for (i=dim; i--; ) {
01077       point= rows[i];
01078       if (point != point0) {
01079         dist= *offset + (point[0]*normal[0] + point[1]*normal[1]
01080                + point[2]*normal[2]);
01081         if (dist > maxround || dist < -maxround) {
01082           *nearzero= True;
01083           break;
01084         }
01085       }
01086     }
01087   }else if (dim == 4) {
01088     normal[0]= - det3_(dY(2,0), dZ(2,0), dW(2,0),
01089                         dY(1,0), dZ(1,0), dW(1,0),
01090                         dY(3,0), dZ(3,0), dW(3,0));
01091     normal[1]=   det3_(dX(2,0), dZ(2,0), dW(2,0),
01092                         dX(1,0), dZ(1,0), dW(1,0),
01093                         dX(3,0), dZ(3,0), dW(3,0));
01094     normal[2]= - det3_(dX(2,0), dY(2,0), dW(2,0),
01095                         dX(1,0), dY(1,0), dW(1,0),
01096                         dX(3,0), dY(3,0), dW(3,0));
01097     normal[3]=   det3_(dX(2,0), dY(2,0), dZ(2,0),
01098                         dX(1,0), dY(1,0), dZ(1,0),
01099                         dX(3,0), dY(3,0), dZ(3,0));
01100     qh_normalize2 (normal, dim, toporient, NULL, NULL);
01101     *offset= -(point0[0]*normal[0] + point0[1]*normal[1]
01102                + point0[2]*normal[2] + point0[3]*normal[3]);
01103     maxround= qh DISTround;
01104     for (i=dim; i--; ) {
01105       point= rows[i];
01106       if (point != point0) {
01107         dist= *offset + (point[0]*normal[0] + point[1]*normal[1]
01108                + point[2]*normal[2] + point[3]*normal[3]);
01109         if (dist > maxround || dist < -maxround) {
01110           *nearzero= True;
01111           break;
01112         }
01113       }
01114     }
01115   }
01116   if (*nearzero) {
01117     zzinc_(Zminnorm);
01118     trace0((qh ferr, "qh_sethyperplane_det: degenerate norm during p%d.\n", qh furthest_id));
01119     zzinc_(Znearlysingular);
01120   }
01121 } /* sethyperplane_det */
01122 
01123 
01124 /*-<a                             href="qh-geom.htm#TOC"
01125   >-------------------------------</a><a name="sethyperplane_gauss">-</a>
01126   
01127   qh_sethyperplane_gauss( dim, rows, point0, toporient, normal, offset, nearzero )
01128     given (dim-1) X dim array of rows[i]= V_{i+1} - V_0 (point0)
01129     set normalized hyperplane equation from oriented simplex
01130 
01131   returns:
01132     normal (normalized)
01133     offset (places point0 on the hyperplane)
01134 
01135   notes:
01136     if nearzero
01137       orientation may be incorrect because of incorrect sign flips in gausselim
01138     solves [V_n-V_0,...,V_1-V_0, 0 .. 0 1] * N == [0 .. 0 1] 
01139         or [V_n-V_0,...,V_1-V_0, 0 .. 0 1] * N == [0] 
01140     i.e., N is normal to the hyperplane, and the unnormalized
01141         distance to [0 .. 1] is either 1 or   0
01142 
01143   design:
01144     perform gaussian elimination
01145     flip sign for negative values
01146     perform back substitution 
01147     normalize result
01148     compute offset
01149 */
01150 void qh_sethyperplane_gauss (int dim, coordT **rows, pointT *point0, 
01151                 boolT toporient, coordT *normal, coordT *offset, boolT *nearzero) {
01152   coordT *pointcoord, *normalcoef;
01153   int k;
01154   boolT sign= toporient, nearzero2= False;
01155   
01156   qh_gausselim(rows, dim-1, dim, &sign, nearzero);
01157   for(k= dim-1; k--; ) {
01158     if ((rows[k])[k] < 0)
01159       sign ^= 1;
01160   }
01161   if (*nearzero) {
01162     zzinc_(Znearlysingular);
01163     trace0((qh ferr, "qh_sethyperplane_gauss: nearly singular or axis parallel hyperplane during p%d.\n", qh furthest_id));
01164     qh_backnormal(rows, dim-1, dim, sign, normal, &nearzero2);
01165   }else {
01166     qh_backnormal(rows, dim-1, dim, sign, normal, &nearzero2);
01167     if (nearzero2) {
01168       zzinc_(Znearlysingular);
01169       trace0((qh ferr, "qh_sethyperplane_gauss: singular or axis parallel hyperplane at normalization during p%d.\n", qh furthest_id));
01170     }
01171   }
01172   if (nearzero2)
01173     *nearzero= True;
01174   qh_normalize2(normal, dim, True, NULL, NULL);
01175   pointcoord= point0;
01176   normalcoef= normal;
01177   *offset= -(*pointcoord++ * *normalcoef++);
01178   for(k= dim-1; k--; )
01179     *offset -= *pointcoord++ * *normalcoef++;
01180 } /* sethyperplane_gauss */
01181 
01182
 
Powered by Plone

This site conforms to the following standards: