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

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