1 /*<html><pre>  -<a                             href="qh-geom_r.htm"
2   >-------------------------------</a><a name="TOP">-</a>
3 
4    geom_r.c
5    geometric routines of qhull
6 
7    see qh-geom_r.htm and geom_r.h
8 
9    Copyright (c) 1993-2019 The Geometry Center.
10    $Id: //main/2019/qhull/src/libqhull_r/geom_r.c#4 $$Change: 2712 $
11    $DateTime: 2019/06/28 12:57:00 $$Author: bbarber $
12 
13    infrequent code goes into geom2_r.c
14 */
15 
16 #include "qhull_ra.h"
17 
18 /*-<a                             href="qh-geom_r.htm#TOC"
19   >-------------------------------</a><a name="distplane">-</a>
20 
21   qh_distplane(qh, point, facet, dist )
22     return distance from point to facet
23 
24   returns:
25     dist
26     if qh.RANDOMdist, joggles result
27 
28   notes:
29     dist > 0 if point is above facet (i.e., outside)
30     does not error (for qh_sortfacets, qh_outerinner)
31     for nearly coplanar points, the returned values may be duplicates
32       for example pairs of nearly incident points, rbox 175 C1,2e-13 t1538759579 | qhull d T4
33       622 qh_distplane: e-014  # count of two or more duplicate values for unique calls
34       258 qh_distplane: e-015
35       38 qh_distplane: e-016
36       40 qh_distplane: e-017
37       6 qh_distplane: e-018
38       5 qh_distplane: -e-018
39       33 qh_distplane: -e-017
40          3153 qh_distplane: -2.775557561562891e-017  # duplicated value for 3153 unique calls
41       42 qh_distplane: -e-016
42       307 qh_distplane: -e-015
43       1271 qh_distplane: -e-014
44       13 qh_distplane: -e-013
45 
46   see:
47     qh_distnorm in geom2_r.c
48     qh_distplane [geom_r.c], QhullFacet::distance, and QhullHyperplane::distance are copies
49 */
qh_distplane(qhT * qh,pointT * point,facetT * facet,realT * dist)50 void qh_distplane(qhT *qh, pointT *point, facetT *facet, realT *dist) {
51   coordT *normal= facet->normal, *coordp, randr;
52   int k;
53 
54   switch (qh->hull_dim){
55   case 2:
56     *dist= facet->offset + point[0] * normal[0] + point[1] * normal[1];
57     break;
58   case 3:
59     *dist= facet->offset + point[0] * normal[0] + point[1] * normal[1] + point[2] * normal[2];
60     break;
61   case 4:
62     *dist= facet->offset+point[0]*normal[0]+point[1]*normal[1]+point[2]*normal[2]+point[3]*normal[3];
63     break;
64   case 5:
65     *dist= facet->offset+point[0]*normal[0]+point[1]*normal[1]+point[2]*normal[2]+point[3]*normal[3]+point[4]*normal[4];
66     break;
67   case 6:
68     *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];
69     break;
70   case 7:
71     *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];
72     break;
73   case 8:
74     *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];
75     break;
76   default:
77     *dist= facet->offset;
78     coordp= point;
79     for (k=qh->hull_dim; k--; )
80       *dist += *coordp++ * *normal++;
81     break;
82   }
83   zzinc_(Zdistplane);
84   if (!qh->RANDOMdist && qh->IStracing < 4)
85     return;
86   if (qh->RANDOMdist) {
87     randr= qh_RANDOMint;
88     *dist += (2.0 * randr / qh_RANDOMmax - 1.0) *
89       qh->RANDOMfactor * qh->MAXabs_coord;
90   }
91 #ifndef qh_NOtrace
92   if (qh->IStracing >= 4) {
93     qh_fprintf(qh, qh->ferr, 8001, "qh_distplane: ");
94     qh_fprintf(qh, qh->ferr, 8002, qh_REAL_1, *dist);
95     qh_fprintf(qh, qh->ferr, 8003, "from p%d to f%d\n", qh_pointid(qh, point), facet->id);
96   }
97 #endif
98   return;
99 } /* distplane */
100 
101 
102 /*-<a                             href="qh-geom_r.htm#TOC"
103   >-------------------------------</a><a name="findbest">-</a>
104 
105   qh_findbest(qh, point, startfacet, bestoutside, qh_ISnewfacets, qh_NOupper, dist, isoutside, numpart )
106     find facet that is furthest below a point
107     for upperDelaunay facets
108       returns facet only if !qh_NOupper and clearly above
109 
110   input:
111     starts search at 'startfacet' (can not be flipped)
112     if !bestoutside(qh_ALL), stops at qh.MINoutside
113 
114   returns:
115     best facet (reports error if NULL)
116     early out if isoutside defined and bestdist > qh.MINoutside
117     dist is distance to facet
118     isoutside is true if point is outside of facet
119     numpart counts the number of distance tests
120 
121   see also:
122     qh_findbestnew()
123 
124   notes:
125     If merging (testhorizon), searches horizon facets of coplanar best facets because
126     after qh_distplane, this and qh_partitionpoint are the most expensive in 3-d
127       avoid calls to distplane, function calls, and real number operations.
128     caller traces result
129     Optimized for outside points.   Tried recording a search set for qh_findhorizon.
130     Made code more complicated.
131 
132   when called by qh_partitionvisible():
133     indicated by qh_ISnewfacets
134     qh.newfacet_list is list of simplicial, new facets
135     qh_findbestnew set if qh_sharpnewfacets returns True (to use qh_findbestnew)
136     qh.bestfacet_notsharp set if qh_sharpnewfacets returns False
137 
138   when called by qh_findfacet(), qh_partitionpoint(), qh_partitioncoplanar(),
139                  qh_check_bestdist(), qh_addpoint()
140     indicated by !qh_ISnewfacets
141     returns best facet in neighborhood of given facet
142       this is best facet overall if dist >= -qh.MAXcoplanar
143         or hull has at least a "spherical" curvature
144 
145   design:
146     initialize and test for early exit
147     repeat while there are better facets
148       for each neighbor of facet
149         exit if outside facet found
150         test for better facet
151     if point is inside and partitioning
152       test for new facets with a "sharp" intersection
153       if so, future calls go to qh_findbestnew()
154     test horizon facets
155 */
qh_findbest(qhT * qh,pointT * point,facetT * startfacet,boolT bestoutside,boolT isnewfacets,boolT noupper,realT * dist,boolT * isoutside,int * numpart)156 facetT *qh_findbest(qhT *qh, pointT *point, facetT *startfacet,
157                      boolT bestoutside, boolT isnewfacets, boolT noupper,
158                      realT *dist, boolT *isoutside, int *numpart) {
159   realT bestdist= -REALmax/2 /* avoid underflow */;
160   facetT *facet, *neighbor, **neighborp;
161   facetT *bestfacet= NULL, *lastfacet= NULL;
162   int oldtrace= qh->IStracing;
163   unsigned int visitid= ++qh->visit_id;
164   int numpartnew=0;
165   boolT testhorizon= True; /* needed if precise, e.g., rbox c D6 | qhull Q0 Tv */
166 
167   zinc_(Zfindbest);
168 #ifndef qh_NOtrace
169   if (qh->IStracing >= 4 || (qh->TRACElevel && qh->TRACEpoint >= 0 && qh->TRACEpoint == qh_pointid(qh, point))) {
170     if (qh->TRACElevel > qh->IStracing)
171       qh->IStracing= qh->TRACElevel;
172     qh_fprintf(qh, qh->ferr, 8004, "qh_findbest: point p%d starting at f%d isnewfacets? %d, unless %d exit if > %2.2g,",
173              qh_pointid(qh, point), startfacet->id, isnewfacets, bestoutside, qh->MINoutside);
174     qh_fprintf(qh, qh->ferr, 8005, " testhorizon? %d, noupper? %d,", testhorizon, noupper);
175     qh_fprintf(qh, qh->ferr, 8006, " Last qh_addpoint p%d,", qh->furthest_id);
176     qh_fprintf(qh, qh->ferr, 8007, " Last merge #%d, max_outside %2.2g\n", zzval_(Ztotmerge), qh->max_outside);
177   }
178 #endif
179   if (isoutside)
180     *isoutside= True;
181   if (!startfacet->flipped) {  /* test startfacet before testing its neighbors */
182     *numpart= 1;
183     qh_distplane(qh, point, startfacet, dist);  /* this code is duplicated below */
184     if (!bestoutside && *dist >= qh->MINoutside
185     && (!startfacet->upperdelaunay || !noupper)) {
186       bestfacet= startfacet;
187       goto LABELreturn_best;
188     }
189     bestdist= *dist;
190     if (!startfacet->upperdelaunay) {
191       bestfacet= startfacet;
192     }
193   }else
194     *numpart= 0;
195   startfacet->visitid= visitid;
196   facet= startfacet;
197   while (facet) {
198     trace4((qh, qh->ferr, 4001, "qh_findbest: neighbors of f%d, bestdist %2.2g f%d\n",
199                 facet->id, bestdist, getid_(bestfacet)));
200     lastfacet= facet;
201     FOREACHneighbor_(facet) {
202       if (!neighbor->newfacet && isnewfacets)
203         continue;
204       if (neighbor->visitid == visitid)
205         continue;
206       neighbor->visitid= visitid;
207       if (!neighbor->flipped) {  /* code duplicated above */
208         (*numpart)++;
209         qh_distplane(qh, point, neighbor, dist);
210         if (*dist > bestdist) {
211           if (!bestoutside && *dist >= qh->MINoutside
212           && (!neighbor->upperdelaunay || !noupper)) {
213             bestfacet= neighbor;
214             goto LABELreturn_best;
215           }
216           if (!neighbor->upperdelaunay) {
217             bestfacet= neighbor;
218             bestdist= *dist;
219             break; /* switch to neighbor */
220           }else if (!bestfacet) {
221             bestdist= *dist;
222             break; /* switch to neighbor */
223           }
224         } /* end of *dist>bestdist */
225       } /* end of !flipped */
226     } /* end of FOREACHneighbor */
227     facet= neighbor;  /* non-NULL only if *dist>bestdist */
228   } /* end of while facet (directed search) */
229   if (isnewfacets) {
230     if (!bestfacet) { /* startfacet is upperdelaunay (or flipped) w/o !flipped newfacet neighbors */
231       bestdist= -REALmax/2;
232       bestfacet= qh_findbestnew(qh, point, qh->newfacet_list, &bestdist, bestoutside, isoutside, &numpartnew);
233       testhorizon= False; /* qh_findbestnew calls qh_findbesthorizon */
234     }else if (!qh->findbest_notsharp && bestdist < -qh->DISTround) {
235       if (qh_sharpnewfacets(qh)) {
236         /* seldom used, qh_findbestnew will retest all facets */
237         zinc_(Zfindnewsharp);
238         bestfacet= qh_findbestnew(qh, point, bestfacet, &bestdist, bestoutside, isoutside, &numpartnew);
239         testhorizon= False; /* qh_findbestnew calls qh_findbesthorizon */
240         qh->findbestnew= True;
241       }else
242         qh->findbest_notsharp= True;
243     }
244   }
245   if (!bestfacet)
246     bestfacet= qh_findbestlower(qh, lastfacet, point, &bestdist, numpart); /* lastfacet is non-NULL because startfacet is non-NULL */
247   if (testhorizon) /* qh_findbestnew not called */
248     bestfacet= qh_findbesthorizon(qh, !qh_IScheckmax, point, bestfacet, noupper, &bestdist, &numpartnew);
249   *dist= bestdist;
250   if (isoutside && bestdist < qh->MINoutside)
251     *isoutside= False;
252 LABELreturn_best:
253   zadd_(Zfindbesttot, *numpart);
254   zmax_(Zfindbestmax, *numpart);
255   (*numpart) += numpartnew;
256   qh->IStracing= oldtrace;
257   return bestfacet;
258 }  /* findbest */
259 
260 
261 /*-<a                             href="qh-geom_r.htm#TOC"
262   >-------------------------------</a><a name="findbesthorizon">-</a>
263 
264   qh_findbesthorizon(qh, qh_IScheckmax, point, startfacet, qh_NOupper, &bestdist, &numpart )
265     search coplanar and better horizon facets from startfacet/bestdist
266     ischeckmax turns off statistics and minsearch update
267     all arguments must be initialized, including *bestdist and *numpart
268     qh.coplanarfacetset used to maintain current search set, reset whenever best facet is substantially better
269   returns(ischeckmax):
270     best facet
271     updates f.maxoutside for neighbors of searched facets (if qh_MAXoutside)
272   returns(!ischeckmax):
273     best facet that is not upperdelaunay or newfacet (qh.first_newfacet)
274     allows upperdelaunay that is clearly outside
275   returns:
276     bestdist is distance to bestfacet
277     numpart -- updates number of distance tests
278 
279   notes:
280     called by qh_findbest if point is not outside a facet (directed search)
281     called by qh_findbestnew if point is not outside a new facet
282     called by qh_check_maxout for each point in hull
283     called by qh_check_bestdist for each point in hull (rarely used)
284 
285     no early out -- use qh_findbest() or qh_findbestnew()
286     Searches coplanar or better horizon facets
287 
288   when called by qh_check_maxout() (qh_IScheckmax)
289     startfacet must be closest to the point
290       Otherwise, if point is beyond and below startfacet, startfacet may be a local minimum
291       even though other facets are below the point.
292     updates facet->maxoutside for good, visited facets
293     may return NULL
294 
295     searchdist is qh.max_outside + 2 * DISTround
296       + max( MINvisible('Vn'), MAXcoplanar('Un'));
297     This setting is a guess.  It must be at least max_outside + 2*DISTround
298     because a facet may have a geometric neighbor across a vertex
299 
300   design:
301     for each horizon facet of coplanar best facets
302       continue if clearly inside
303       unless upperdelaunay or clearly outside
304          update best facet
305 */
qh_findbesthorizon(qhT * qh,boolT ischeckmax,pointT * point,facetT * startfacet,boolT noupper,realT * bestdist,int * numpart)306 facetT *qh_findbesthorizon(qhT *qh, boolT ischeckmax, pointT* point, facetT *startfacet, boolT noupper, realT *bestdist, int *numpart) {
307   facetT *bestfacet= startfacet;
308   realT dist;
309   facetT *neighbor, **neighborp, *facet;
310   facetT *nextfacet= NULL; /* optimize last facet of coplanarfacetset */
311   int numpartinit= *numpart, coplanarfacetset_size, numcoplanar= 0, numfacet= 0;
312   unsigned int visitid= ++qh->visit_id;
313   boolT newbest= False; /* for tracing */
314   realT minsearch, searchdist;  /* skip facets that are too far from point */
315   boolT is_5x_minsearch;
316 
317   if (!ischeckmax) {
318     zinc_(Zfindhorizon);
319   }else {
320 #if qh_MAXoutside
321     if ((!qh->ONLYgood || startfacet->good) && *bestdist > startfacet->maxoutside)
322       startfacet->maxoutside= *bestdist;
323 #endif
324   }
325   searchdist= qh_SEARCHdist; /* an expression, a multiple of qh.max_outside and precision constants */
326   minsearch= *bestdist - searchdist;
327   if (ischeckmax) {
328     /* Always check coplanar facets.  Needed for RBOX 1000 s Z1 G1e-13 t996564279 | QHULL Tv */
329     minimize_(minsearch, -searchdist);
330   }
331   coplanarfacetset_size= 0;
332   startfacet->visitid= visitid;
333   facet= startfacet;
334   while (True) {
335     numfacet++;
336     is_5x_minsearch= (ischeckmax && facet->nummerge > 10 && qh_setsize(qh, facet->neighbors) > 100);  /* QH11033 FIX: qh_findbesthorizon: many tests for facets with many merges and neighbors.  Can hide coplanar facets, e.g., 'rbox 1000 s Z1 G1e-13' with 4400+ neighbors */
337     trace4((qh, qh->ferr, 4002, "qh_findbesthorizon: test neighbors of f%d bestdist %2.2g f%d ischeckmax? %d noupper? %d minsearch %2.2g is_5x? %d searchdist %2.2g\n",
338                 facet->id, *bestdist, getid_(bestfacet), ischeckmax, noupper,
339                 minsearch, is_5x_minsearch, searchdist));
340     FOREACHneighbor_(facet) {
341       if (neighbor->visitid == visitid)
342         continue;
343       neighbor->visitid= visitid;
344       if (!neighbor->flipped) {  /* neighbors of flipped facets always searched via nextfacet */
345         qh_distplane(qh, point, neighbor, &dist); /* duplicate qh_distpane for new facets, they may be coplanar */
346         (*numpart)++;
347         if (dist > *bestdist) {
348           if (!neighbor->upperdelaunay || ischeckmax || (!noupper && dist >= qh->MINoutside)) {
349             if (!ischeckmax) {
350               minsearch= dist - searchdist;
351               if (dist > *bestdist + searchdist) {
352                 zinc_(Zfindjump);  /* everything in qh.coplanarfacetset at least searchdist below */
353                 coplanarfacetset_size= 0;
354               }
355             }
356             bestfacet= neighbor;
357             *bestdist= dist;
358             newbest= True;
359           }
360         }else if (is_5x_minsearch) {
361           if (dist < 5 * minsearch)
362             continue; /* skip this neighbor, do not set nextfacet.  dist is negative */
363         }else if (dist < minsearch)
364           continue;  /* skip this neighbor, do not set nextfacet.  If ischeckmax, dist can't be positive */
365 #if qh_MAXoutside
366         if (ischeckmax && dist > neighbor->maxoutside)
367           neighbor->maxoutside= dist;
368 #endif
369       } /* end of !flipped, need to search neighbor */
370       if (nextfacet) {
371         numcoplanar++;
372         if (!coplanarfacetset_size++) {
373           SETfirst_(qh->coplanarfacetset)= nextfacet;
374           SETtruncate_(qh->coplanarfacetset, 1);
375         }else
376           qh_setappend(qh, &qh->coplanarfacetset, nextfacet); /* Was needed for RBOX 1000 s W1e-13 P0 t996547055 | QHULL d Qbb Qc Tv
377                                                  and RBOX 1000 s Z1 G1e-13 t996564279 | qhull Tv  */
378       }
379       nextfacet= neighbor;
380     } /* end of EACHneighbor */
381     facet= nextfacet;
382     if (facet)
383       nextfacet= NULL;
384     else if (!coplanarfacetset_size)
385       break;
386     else if (!--coplanarfacetset_size) {
387       facet= SETfirstt_(qh->coplanarfacetset, facetT);
388       SETtruncate_(qh->coplanarfacetset, 0);
389     }else
390       facet= (facetT *)qh_setdellast(qh->coplanarfacetset);
391   } /* while True, i.e., "for each facet in qh.coplanarfacetset" */
392   if (!ischeckmax) {
393     zadd_(Zfindhorizontot, *numpart - numpartinit);
394     zmax_(Zfindhorizonmax, *numpart - numpartinit);
395     if (newbest)
396       zinc_(Znewbesthorizon);
397   }
398   trace4((qh, qh->ferr, 4003, "qh_findbesthorizon: p%d, newbest? %d, bestfacet f%d, bestdist %2.2g, numfacet %d, coplanarfacets %d, numdist %d\n",
399     qh_pointid(qh, point), newbest, getid_(bestfacet), *bestdist, numfacet, numcoplanar, *numpart - numpartinit));
400   return bestfacet;
401 }  /* findbesthorizon */
402 
403 /*-<a                             href="qh-geom_r.htm#TOC"
404   >-------------------------------</a><a name="findbestnew">-</a>
405 
406   qh_findbestnew(qh, point, startfacet, dist, isoutside, numpart )
407     find best newfacet for point
408     searches all of qh.newfacet_list starting at startfacet
409     searches horizon facets of coplanar best newfacets
410     searches all facets if startfacet == qh.facet_list
411   returns:
412     best new or horizon facet that is not upperdelaunay
413     early out if isoutside and not 'Qf'
414     dist is distance to facet
415     isoutside is true if point is outside of facet
416     numpart is number of distance tests
417 
418   notes:
419     Always used for merged new facets (see qh_USEfindbestnew)
420     Avoids upperdelaunay facet unless (isoutside and outside)
421 
422     Uses qh.visit_id, qh.coplanarfacetset.
423     If share visit_id with qh_findbest, coplanarfacetset is incorrect.
424 
425     If merging (testhorizon), searches horizon facets of coplanar best facets because
426     a point maybe coplanar to the bestfacet, below its horizon facet,
427     and above a horizon facet of a coplanar newfacet.  For example,
428       rbox 1000 s Z1 G1e-13 | qhull
429       rbox 1000 s W1e-13 P0 t992110337 | QHULL d Qbb Qc
430 
431     qh_findbestnew() used if
432        qh_sharpnewfacets -- newfacets contains a sharp angle
433        if many merges, qh_premerge found a merge, or 'Qf' (qh.findbestnew)
434 
435   see also:
436     qh_partitionall() and qh_findbest()
437 
438   design:
439     for each new facet starting from startfacet
440       test distance from point to facet
441       return facet if clearly outside
442       unless upperdelaunay and a lowerdelaunay exists
443          update best facet
444     test horizon facets
445 */
qh_findbestnew(qhT * qh,pointT * point,facetT * startfacet,realT * dist,boolT bestoutside,boolT * isoutside,int * numpart)446 facetT *qh_findbestnew(qhT *qh, pointT *point, facetT *startfacet,
447            realT *dist, boolT bestoutside, boolT *isoutside, int *numpart) {
448   realT bestdist= -REALmax/2;
449   facetT *bestfacet= NULL, *facet;
450   int oldtrace= qh->IStracing, i;
451   unsigned int visitid= ++qh->visit_id;
452   realT distoutside= 0.0;
453   boolT isdistoutside; /* True if distoutside is defined */
454   boolT testhorizon= True; /* needed if precise, e.g., rbox c D6 | qhull Q0 Tv */
455 
456   if (!startfacet || !startfacet->next) {
457     if (qh->MERGING) {
458       qh_fprintf(qh, qh->ferr, 6001, "qhull topology error (qh_findbestnew): merging has formed and deleted a cone of new facets.  Can not continue.\n");
459       qh_errexit(qh, qh_ERRtopology, NULL, NULL);
460     }else {
461       qh_fprintf(qh, qh->ferr, 6002, "qhull internal error (qh_findbestnew): no new facets for point p%d\n",
462               qh->furthest_id);
463       qh_errexit(qh, qh_ERRqhull, NULL, NULL);
464     }
465   }
466   zinc_(Zfindnew);
467   if (qh->BESToutside || bestoutside)
468     isdistoutside= False;
469   else {
470     isdistoutside= True;
471     distoutside= qh_DISToutside; /* multiple of qh.MINoutside & qh.max_outside, see user_r.h */
472   }
473   if (isoutside)
474     *isoutside= True;
475   *numpart= 0;
476 #ifndef qh_NOtrace
477   if (qh->IStracing >= 4 || (qh->TRACElevel && qh->TRACEpoint >= 0 && qh->TRACEpoint == qh_pointid(qh, point))) {
478     if (qh->TRACElevel > qh->IStracing)
479       qh->IStracing= qh->TRACElevel;
480     qh_fprintf(qh, qh->ferr, 8008, "qh_findbestnew: point p%d facet f%d. Stop? %d if dist > %2.2g,",
481              qh_pointid(qh, point), startfacet->id, isdistoutside, distoutside);
482     qh_fprintf(qh, qh->ferr, 8009, " Last qh_addpoint p%d, qh.visit_id %d, vertex_visit %d,",  qh->furthest_id, visitid, qh->vertex_visit);
483     qh_fprintf(qh, qh->ferr, 8010, " Last merge #%d\n", zzval_(Ztotmerge));
484   }
485 #endif
486   /* visit all new facets starting with startfacet, maybe qh->facet_list */
487   for (i=0, facet=startfacet; i < 2; i++, facet= qh->newfacet_list) {
488     FORALLfacet_(facet) {
489       if (facet == startfacet && i)
490         break;
491       facet->visitid= visitid;
492       if (!facet->flipped) {
493         qh_distplane(qh, point, facet, dist);
494         (*numpart)++;
495         if (*dist > bestdist) {
496           if (!facet->upperdelaunay || *dist >= qh->MINoutside) {
497             bestfacet= facet;
498             if (isdistoutside && *dist >= distoutside)
499               goto LABELreturn_bestnew;
500             bestdist= *dist;
501           }
502         }
503       } /* end of !flipped */
504     } /* FORALLfacet from startfacet or qh->newfacet_list */
505   }
506   if (testhorizon || !bestfacet) /* testhorizon is always True.  Keep the same code as qh_findbest */
507     bestfacet= qh_findbesthorizon(qh, !qh_IScheckmax, point, bestfacet ? bestfacet : startfacet,
508                                         !qh_NOupper, &bestdist, numpart);
509   *dist= bestdist;
510   if (isoutside && *dist < qh->MINoutside)
511     *isoutside= False;
512 LABELreturn_bestnew:
513   zadd_(Zfindnewtot, *numpart);
514   zmax_(Zfindnewmax, *numpart);
515   trace4((qh, qh->ferr, 4004, "qh_findbestnew: bestfacet f%d bestdist %2.2g for p%d f%d bestoutside? %d \n",
516     getid_(bestfacet), *dist, qh_pointid(qh, point), startfacet->id, bestoutside));
517   qh->IStracing= oldtrace;
518   return bestfacet;
519 }  /* findbestnew */
520 
521 /* ============ hyperplane functions -- keep code together [?] ============ */
522 
523 /*-<a                             href="qh-geom_r.htm#TOC"
524   >-------------------------------</a><a name="backnormal">-</a>
525 
526   qh_backnormal(qh, rows, numrow, numcol, sign, normal, nearzero )
527     given an upper-triangular rows array and a sign,
528     solve for normal equation x using back substitution over rows U
529 
530   returns:
531      normal= x
532 
533      if will not be able to divzero() when normalized(qh.MINdenom_2 and qh.MINdenom_1_2),
534        if fails on last row
535          this means that the hyperplane intersects [0,..,1]
536          sets last coordinate of normal to sign
537        otherwise
538          sets tail of normal to [...,sign,0,...], i.e., solves for b= [0...0]
539          sets nearzero
540 
541   notes:
542      assumes numrow == numcol-1
543 
544      see Golub & van Loan, 1983, Eq. 4.4-9 for "Gaussian elimination with complete pivoting"
545 
546      solves Ux=b where Ax=b and PA=LU
547      b= [0,...,0,sign or 0]  (sign is either -1 or +1)
548      last row of A= [0,...,0,1]
549 
550      1) Ly=Pb == y=b since P only permutes the 0's of   b
551 
552   design:
553     for each row from end
554       perform back substitution
555       if near zero
556         use qh_divzero for division
557         if zero divide and not last row
558           set tail of normal to 0
559 */
qh_backnormal(qhT * qh,realT ** rows,int numrow,int numcol,boolT sign,coordT * normal,boolT * nearzero)560 void qh_backnormal(qhT *qh, realT **rows, int numrow, int numcol, boolT sign,
561         coordT *normal, boolT *nearzero) {
562   int i, j;
563   coordT *normalp, *normal_tail, *ai, *ak;
564   realT diagonal;
565   boolT waszero;
566   int zerocol= -1;
567 
568   normalp= normal + numcol - 1;
569   *normalp--= (sign ? -1.0 : 1.0);
570   for (i=numrow; i--; ) {
571     *normalp= 0.0;
572     ai= rows[i] + i + 1;
573     ak= normalp+1;
574     for (j=i+1; j < numcol; j++)
575       *normalp -= *ai++ * *ak++;
576     diagonal= (rows[i])[i];
577     if (fabs_(diagonal) > qh->MINdenom_2)
578       *(normalp--) /= diagonal;
579     else {
580       waszero= False;
581       *normalp= qh_divzero(*normalp, diagonal, qh->MINdenom_1_2, &waszero);
582       if (waszero) {
583         zerocol= i;
584         *(normalp--)= (sign ? -1.0 : 1.0);
585         for (normal_tail= normalp+2; normal_tail < normal + numcol; normal_tail++)
586           *normal_tail= 0.0;
587       }else
588         normalp--;
589     }
590   }
591   if (zerocol != -1) {
592     *nearzero= True;
593     trace4((qh, qh->ferr, 4005, "qh_backnormal: zero diagonal at column %d.\n", i));
594     zzinc_(Zback0);
595     qh_joggle_restart(qh, "zero diagonal on back substitution");
596   }
597 } /* backnormal */
598 
599 /*-<a                             href="qh-geom_r.htm#TOC"
600   >-------------------------------</a><a name="gausselim">-</a>
601 
602   qh_gausselim(qh, rows, numrow, numcol, sign )
603     Gaussian elimination with partial pivoting
604 
605   returns:
606     rows is upper triangular (includes row exchanges)
607     flips sign for each row exchange
608     sets nearzero if pivot[k] < qh.NEARzero[k], else clears it
609 
610   notes:
611     if nearzero, the determinant's sign may be incorrect.
612     assumes numrow <= numcol
613 
614   design:
615     for each row
616       determine pivot and exchange rows if necessary
617       test for near zero
618       perform gaussian elimination step
619 */
qh_gausselim(qhT * qh,realT ** rows,int numrow,int numcol,boolT * sign,boolT * nearzero)620 void qh_gausselim(qhT *qh, realT **rows, int numrow, int numcol, boolT *sign, boolT *nearzero) {
621   realT *ai, *ak, *rowp, *pivotrow;
622   realT n, pivot, pivot_abs= 0.0, temp;
623   int i, j, k, pivoti, flip=0;
624 
625   *nearzero= False;
626   for (k=0; k < numrow; k++) {
627     pivot_abs= fabs_((rows[k])[k]);
628     pivoti= k;
629     for (i=k+1; i < numrow; i++) {
630       if ((temp= fabs_((rows[i])[k])) > pivot_abs) {
631         pivot_abs= temp;
632         pivoti= i;
633       }
634     }
635     if (pivoti != k) {
636       rowp= rows[pivoti];
637       rows[pivoti]= rows[k];
638       rows[k]= rowp;
639       *sign ^= 1;
640       flip ^= 1;
641     }
642     if (pivot_abs <= qh->NEARzero[k]) {
643       *nearzero= True;
644       if (pivot_abs == 0.0) {   /* remainder of column == 0 */
645 #ifndef qh_NOtrace
646         if (qh->IStracing >= 4) {
647           qh_fprintf(qh, qh->ferr, 8011, "qh_gausselim: 0 pivot at column %d. (%2.2g < %2.2g)\n", k, pivot_abs, qh->DISTround);
648           qh_printmatrix(qh, qh->ferr, "Matrix:", rows, numrow, numcol);
649         }
650 #endif
651         zzinc_(Zgauss0);
652         qh_joggle_restart(qh, "zero pivot for Gaussian elimination");
653         goto LABELnextcol;
654       }
655     }
656     pivotrow= rows[k] + k;
657     pivot= *pivotrow++;  /* signed value of pivot, and remainder of row */
658     for (i=k+1; i < numrow; i++) {
659       ai= rows[i] + k;
660       ak= pivotrow;
661       n= (*ai++)/pivot;   /* divzero() not needed since |pivot| >= |*ai| */
662       for (j= numcol - (k+1); j--; )
663         *ai++ -= n * *ak++;
664     }
665   LABELnextcol:
666     ;
667   }
668   wmin_(Wmindenom, pivot_abs);  /* last pivot element */
669   if (qh->IStracing >= 5)
670     qh_printmatrix(qh, qh->ferr, "qh_gausselem: result", rows, numrow, numcol);
671 } /* gausselim */
672 
673 
674 /*-<a                             href="qh-geom_r.htm#TOC"
675   >-------------------------------</a><a name="getangle">-</a>
676 
677   qh_getangle(qh, vect1, vect2 )
678     returns the dot product of two vectors
679     if qh.RANDOMdist, joggles result
680 
681   notes:
682     the angle may be > 1.0 or < -1.0 because of roundoff errors
683 
684 */
qh_getangle(qhT * qh,pointT * vect1,pointT * vect2)685 realT qh_getangle(qhT *qh, pointT *vect1, pointT *vect2) {
686   realT angle= 0, randr;
687   int k;
688 
689   for (k=qh->hull_dim; k--; )
690     angle += *vect1++ * *vect2++;
691   if (qh->RANDOMdist) {
692     randr= qh_RANDOMint;
693     angle += (2.0 * randr / qh_RANDOMmax - 1.0) *
694       qh->RANDOMfactor;
695   }
696   trace4((qh, qh->ferr, 4006, "qh_getangle: %4.4g\n", angle));
697   return(angle);
698 } /* getangle */
699 
700 
701 /*-<a                             href="qh-geom_r.htm#TOC"
702   >-------------------------------</a><a name="getcenter">-</a>
703 
704   qh_getcenter(qh, vertices )
705     returns arithmetic center of a set of vertices as a new point
706 
707   notes:
708     allocates point array for center
709 */
qh_getcenter(qhT * qh,setT * vertices)710 pointT *qh_getcenter(qhT *qh, setT *vertices) {
711   int k;
712   pointT *center, *coord;
713   vertexT *vertex, **vertexp;
714   int count= qh_setsize(qh, vertices);
715 
716   if (count < 2) {
717     qh_fprintf(qh, qh->ferr, 6003, "qhull internal error (qh_getcenter): not defined for %d points\n", count);
718     qh_errexit(qh, qh_ERRqhull, NULL, NULL);
719   }
720   center= (pointT *)qh_memalloc(qh, qh->normal_size);
721   for (k=0; k < qh->hull_dim; k++) {
722     coord= center+k;
723     *coord= 0.0;
724     FOREACHvertex_(vertices)
725       *coord += vertex->point[k];
726     *coord /= count;  /* count>=2 by QH6003 */
727   }
728   return(center);
729 } /* getcenter */
730 
731 
732 /*-<a                             href="qh-geom_r.htm#TOC"
733   >-------------------------------</a><a name="getcentrum">-</a>
734 
735   qh_getcentrum(qh, facet )
736     returns the centrum for a facet as a new point
737 
738   notes:
739     allocates the centrum
740 */
qh_getcentrum(qhT * qh,facetT * facet)741 pointT *qh_getcentrum(qhT *qh, facetT *facet) {
742   realT dist;
743   pointT *centrum, *point;
744 
745   point= qh_getcenter(qh, facet->vertices);
746   zzinc_(Zcentrumtests);
747   qh_distplane(qh, point, facet, &dist);
748   centrum= qh_projectpoint(qh, point, facet, dist);
749   qh_memfree(qh, point, qh->normal_size);
750   trace4((qh, qh->ferr, 4007, "qh_getcentrum: for f%d, %d vertices dist= %2.2g\n",
751           facet->id, qh_setsize(qh, facet->vertices), dist));
752   return centrum;
753 } /* getcentrum */
754 
755 
756 /*-<a                             href="qh-geom_r.htm#TOC"
757   >-------------------------------</a><a name="getdistance">-</a>
758 
759   qh_getdistance(qh, facet, neighbor, mindist, maxdist )
760     returns the min and max distance to neighbor of non-neighbor vertices in facet
761 
762   returns:
763     the max absolute value
764 
765   design:
766     for each vertex of facet that is not in neighbor
767       test the distance from vertex to neighbor
768 */
qh_getdistance(qhT * qh,facetT * facet,facetT * neighbor,coordT * mindist,coordT * maxdist)769 coordT qh_getdistance(qhT *qh, facetT *facet, facetT *neighbor, coordT *mindist, coordT *maxdist) {
770   vertexT *vertex, **vertexp;
771   coordT dist, maxd, mind;
772 
773   FOREACHvertex_(facet->vertices)
774     vertex->seen= False;
775   FOREACHvertex_(neighbor->vertices)
776     vertex->seen= True;
777   mind= 0.0;
778   maxd= 0.0;
779   FOREACHvertex_(facet->vertices) {
780     if (!vertex->seen) {
781       zzinc_(Zbestdist);
782       qh_distplane(qh, vertex->point, neighbor, &dist);
783       if (dist < mind)
784         mind= dist;
785       else if (dist > maxd)
786         maxd= dist;
787     }
788   }
789   *mindist= mind;
790   *maxdist= maxd;
791   mind= -mind;
792   if (maxd > mind)
793     return maxd;
794   else
795     return mind;
796 } /* getdistance */
797 
798 
799 /*-<a                             href="qh-geom_r.htm#TOC"
800   >-------------------------------</a><a name="normalize">-</a>
801 
802   qh_normalize(qh, normal, dim, toporient )
803     normalize a vector and report if too small
804     does not use min norm
805 
806   see:
807     qh_normalize2
808 */
qh_normalize(qhT * qh,coordT * normal,int dim,boolT toporient)809 void qh_normalize(qhT *qh, coordT *normal, int dim, boolT toporient) {
810   qh_normalize2(qh, normal, dim, toporient, NULL, NULL);
811 } /* normalize */
812 
813 /*-<a                             href="qh-geom_r.htm#TOC"
814   >-------------------------------</a><a name="normalize2">-</a>
815 
816   qh_normalize2(qh, normal, dim, toporient, minnorm, ismin )
817     normalize a vector and report if too small
818     qh.MINdenom/MINdenom1 are the upper limits for divide overflow
819 
820   returns:
821     normalized vector
822     flips sign if !toporient
823     if minnorm non-NULL,
824       sets ismin if normal < minnorm
825 
826   notes:
827     if zero norm
828        sets all elements to sqrt(1.0/dim)
829     if divide by zero (divzero())
830        sets largest element to   +/-1
831        bumps Znearlysingular
832 
833   design:
834     computes norm
835     test for minnorm
836     if not near zero
837       normalizes normal
838     else if zero norm
839       sets normal to standard value
840     else
841       uses qh_divzero to normalize
842       if nearzero
843         sets norm to direction of maximum value
844 */
qh_normalize2(qhT * qh,coordT * normal,int dim,boolT toporient,realT * minnorm,boolT * ismin)845 void qh_normalize2(qhT *qh, coordT *normal, int dim, boolT toporient,
846             realT *minnorm, boolT *ismin) {
847   int k;
848   realT *colp, *maxp, norm= 0, temp, *norm1, *norm2, *norm3;
849   boolT zerodiv;
850 
851   norm1= normal+1;
852   norm2= normal+2;
853   norm3= normal+3;
854   if (dim == 2)
855     norm= sqrt((*normal)*(*normal) + (*norm1)*(*norm1));
856   else if (dim == 3)
857     norm= sqrt((*normal)*(*normal) + (*norm1)*(*norm1) + (*norm2)*(*norm2));
858   else if (dim == 4) {
859     norm= sqrt((*normal)*(*normal) + (*norm1)*(*norm1) + (*norm2)*(*norm2)
860                + (*norm3)*(*norm3));
861   }else if (dim > 4) {
862     norm= (*normal)*(*normal) + (*norm1)*(*norm1) + (*norm2)*(*norm2)
863                + (*norm3)*(*norm3);
864     for (k=dim-4, colp=normal+4; k--; colp++)
865       norm += (*colp) * (*colp);
866     norm= sqrt(norm);
867   }
868   if (minnorm) {
869     if (norm < *minnorm)
870       *ismin= True;
871     else
872       *ismin= False;
873   }
874   wmin_(Wmindenom, norm);
875   if (norm > qh->MINdenom) {
876     if (!toporient)
877       norm= -norm;
878     *normal /= norm;
879     *norm1 /= norm;
880     if (dim == 2)
881       ; /* all done */
882     else if (dim == 3)
883       *norm2 /= norm;
884     else if (dim == 4) {
885       *norm2 /= norm;
886       *norm3 /= norm;
887     }else if (dim >4) {
888       *norm2 /= norm;
889       *norm3 /= norm;
890       for (k=dim-4, colp=normal+4; k--; )
891         *colp++ /= norm;
892     }
893   }else if (norm == 0.0) {
894     temp= sqrt(1.0/dim);
895     for (k=dim, colp=normal; k--; )
896       *colp++= temp;
897   }else {
898     if (!toporient)
899       norm= -norm;
900     for (k=dim, colp=normal; k--; colp++) { /* k used below */
901       temp= qh_divzero(*colp, norm, qh->MINdenom_1, &zerodiv);
902       if (!zerodiv)
903         *colp= temp;
904       else {
905         maxp= qh_maxabsval(normal, dim);
906         temp= ((*maxp * norm >= 0.0) ? 1.0 : -1.0);
907         for (k=dim, colp=normal; k--; colp++)
908           *colp= 0.0;
909         *maxp= temp;
910         zzinc_(Znearlysingular);
911         /* qh_joggle_restart ignored for Znearlysingular, normal part of qh_sethyperplane_gauss */
912         trace0((qh, qh->ferr, 1, "qh_normalize: norm=%2.2g too small during p%d\n",
913                norm, qh->furthest_id));
914         return;
915       }
916     }
917   }
918 } /* normalize */
919 
920 
921 /*-<a                             href="qh-geom_r.htm#TOC"
922   >-------------------------------</a><a name="projectpoint">-</a>
923 
924   qh_projectpoint(qh, point, facet, dist )
925     project point onto a facet by dist
926 
927   returns:
928     returns a new point
929 
930   notes:
931     if dist= distplane(point,facet)
932       this projects point to hyperplane
933     assumes qh_memfree_() is valid for normal_size
934 */
qh_projectpoint(qhT * qh,pointT * point,facetT * facet,realT dist)935 pointT *qh_projectpoint(qhT *qh, pointT *point, facetT *facet, realT dist) {
936   pointT *newpoint, *np, *normal;
937   int normsize= qh->normal_size;
938   int k;
939   void **freelistp; /* used if !qh_NOmem by qh_memalloc_() */
940 
941   qh_memalloc_(qh, normsize, freelistp, newpoint, pointT);
942   np= newpoint;
943   normal= facet->normal;
944   for (k=qh->hull_dim; k--; )
945     *(np++)= *point++ - dist * *normal++;
946   return(newpoint);
947 } /* projectpoint */
948 
949 
950 /*-<a                             href="qh-geom_r.htm#TOC"
951   >-------------------------------</a><a name="setfacetplane">-</a>
952 
953   qh_setfacetplane(qh, facet )
954     sets the hyperplane for a facet
955     if qh.RANDOMdist, joggles hyperplane
956 
957   notes:
958     uses global buffers qh.gm_matrix and qh.gm_row
959     overwrites facet->normal if already defined
960     updates Wnewvertex if PRINTstatistics
961     sets facet->upperdelaunay if upper envelope of Delaunay triangulation
962 
963   design:
964     copy vertex coordinates to qh.gm_matrix/gm_row
965     compute determinate
966     if nearzero
967       recompute determinate with gaussian elimination
968       if nearzero
969         force outside orientation by testing interior point
970 */
qh_setfacetplane(qhT * qh,facetT * facet)971 void qh_setfacetplane(qhT *qh, facetT *facet) {
972   pointT *point;
973   vertexT *vertex, **vertexp;
974   int normsize= qh->normal_size;
975   int k,i, oldtrace= 0;
976   realT dist;
977   void **freelistp; /* used if !qh_NOmem by qh_memalloc_() */
978   coordT *coord, *gmcoord;
979   pointT *point0= SETfirstt_(facet->vertices, vertexT)->point;
980   boolT nearzero= False;
981 
982   zzinc_(Zsetplane);
983   if (!facet->normal)
984     qh_memalloc_(qh, normsize, freelistp, facet->normal, coordT);
985 #ifndef qh_NOtrace
986   if (facet == qh->tracefacet) {
987     oldtrace= qh->IStracing;
988     qh->IStracing= 5;
989     qh_fprintf(qh, qh->ferr, 8012, "qh_setfacetplane: facet f%d created.\n", facet->id);
990     qh_fprintf(qh, qh->ferr, 8013, "  Last point added to hull was p%d.", qh->furthest_id);
991     if (zzval_(Ztotmerge))
992       qh_fprintf(qh, qh->ferr, 8014, "  Last merge was #%d.", zzval_(Ztotmerge));
993     qh_fprintf(qh, qh->ferr, 8015, "\n\nCurrent summary is:\n");
994       qh_printsummary(qh, qh->ferr);
995   }
996 #endif
997   if (qh->hull_dim <= 4) {
998     i= 0;
999     if (qh->RANDOMdist) {
1000       gmcoord= qh->gm_matrix;
1001       FOREACHvertex_(facet->vertices) {
1002         qh->gm_row[i++]= gmcoord;
1003         coord= vertex->point;
1004         for (k=qh->hull_dim; k--; )
1005           *(gmcoord++)= *coord++ * qh_randomfactor(qh, qh->RANDOMa, qh->RANDOMb);
1006       }
1007     }else {
1008       FOREACHvertex_(facet->vertices)
1009        qh->gm_row[i++]= vertex->point;
1010     }
1011     qh_sethyperplane_det(qh, qh->hull_dim, qh->gm_row, point0, facet->toporient,
1012                 facet->normal, &facet->offset, &nearzero);
1013   }
1014   if (qh->hull_dim > 4 || nearzero) {
1015     i= 0;
1016     gmcoord= qh->gm_matrix;
1017     FOREACHvertex_(facet->vertices) {
1018       if (vertex->point != point0) {
1019         qh->gm_row[i++]= gmcoord;
1020         coord= vertex->point;
1021         point= point0;
1022         for (k=qh->hull_dim; k--; )
1023           *(gmcoord++)= *coord++ - *point++;
1024       }
1025     }
1026     qh->gm_row[i]= gmcoord;  /* for areasimplex */
1027     if (qh->RANDOMdist) {
1028       gmcoord= qh->gm_matrix;
1029       for (i=qh->hull_dim-1; i--; ) {
1030         for (k=qh->hull_dim; k--; )
1031           *(gmcoord++) *= qh_randomfactor(qh, qh->RANDOMa, qh->RANDOMb);
1032       }
1033     }
1034     qh_sethyperplane_gauss(qh, qh->hull_dim, qh->gm_row, point0, facet->toporient,
1035                 facet->normal, &facet->offset, &nearzero);
1036     if (nearzero) {
1037       if (qh_orientoutside(qh, facet)) {
1038         trace0((qh, qh->ferr, 2, "qh_setfacetplane: flipped orientation due to nearzero gauss and interior_point test.  During p%d\n", qh->furthest_id));
1039       /* this is part of using Gaussian Elimination.  For example in 5-d
1040            1 1 1 1 0
1041            1 1 1 1 1
1042            0 0 0 1 0
1043            0 1 0 0 0
1044            1 0 0 0 0
1045            norm= 0.38 0.38 -0.76 0.38 0
1046          has a determinate of 1, but g.e. after subtracting pt. 0 has
1047          0's in the diagonal, even with full pivoting.  It does work
1048          if you subtract pt. 4 instead. */
1049       }
1050     }
1051   }
1052   facet->upperdelaunay= False;
1053   if (qh->DELAUNAY) {
1054     if (qh->UPPERdelaunay) {     /* matches qh_triangulate_facet and qh.lower_threshold in qh_initbuild */
1055       if (facet->normal[qh->hull_dim -1] >= qh->ANGLEround * qh_ZEROdelaunay)
1056         facet->upperdelaunay= True;
1057     }else {
1058       if (facet->normal[qh->hull_dim -1] > -qh->ANGLEround * qh_ZEROdelaunay)
1059         facet->upperdelaunay= True;
1060     }
1061   }
1062   if (qh->PRINTstatistics || qh->IStracing || qh->TRACElevel || qh->JOGGLEmax < REALmax) {
1063     qh->old_randomdist= qh->RANDOMdist;
1064     qh->RANDOMdist= False;
1065     FOREACHvertex_(facet->vertices) {
1066       if (vertex->point != point0) {
1067         boolT istrace= False;
1068         zinc_(Zdiststat);
1069         qh_distplane(qh, vertex->point, facet, &dist);
1070         dist= fabs_(dist);
1071         zinc_(Znewvertex);
1072         wadd_(Wnewvertex, dist);
1073         if (dist > wwval_(Wnewvertexmax)) {
1074           wwval_(Wnewvertexmax)= dist;
1075           if (dist > qh->max_outside) {
1076             qh->max_outside= dist;  /* used by qh_maxouter(qh) */
1077             if (dist > qh->TRACEdist)
1078               istrace= True;
1079           }
1080         }else if (-dist > qh->TRACEdist)
1081           istrace= True;
1082         if (istrace) {
1083           qh_fprintf(qh, qh->ferr, 3060, "qh_setfacetplane: ====== vertex p%d(v%d) increases max_outside to %2.2g for new facet f%d last p%d\n",
1084                 qh_pointid(qh, vertex->point), vertex->id, dist, facet->id, qh->furthest_id);
1085           qh_errprint(qh, "DISTANT", facet, NULL, NULL, NULL);
1086         }
1087       }
1088     }
1089     qh->RANDOMdist= qh->old_randomdist;
1090   }
1091 #ifndef qh_NOtrace
1092   if (qh->IStracing >= 4) {
1093     qh_fprintf(qh, qh->ferr, 8017, "qh_setfacetplane: f%d offset %2.2g normal: ",
1094              facet->id, facet->offset);
1095     for (k=0; k < qh->hull_dim; k++)
1096       qh_fprintf(qh, qh->ferr, 8018, "%2.2g ", facet->normal[k]);
1097     qh_fprintf(qh, qh->ferr, 8019, "\n");
1098   }
1099 #endif
1100   qh_checkflipped(qh, facet, NULL, qh_ALL);
1101   if (facet == qh->tracefacet) {
1102     qh->IStracing= oldtrace;
1103     qh_printfacet(qh, qh->ferr, facet);
1104   }
1105 } /* setfacetplane */
1106 
1107 
1108 /*-<a                             href="qh-geom_r.htm#TOC"
1109   >-------------------------------</a><a name="sethyperplane_det">-</a>
1110 
1111   qh_sethyperplane_det(qh, dim, rows, point0, toporient, normal, offset, nearzero )
1112     given dim X dim array indexed by rows[], one row per point,
1113         toporient(flips all signs),
1114         and point0 (any row)
1115     set normalized hyperplane equation from oriented simplex
1116 
1117   returns:
1118     normal (normalized)
1119     offset (places point0 on the hyperplane)
1120     sets nearzero if hyperplane not through points
1121 
1122   notes:
1123     only defined for dim == 2..4
1124     rows[] is not modified
1125     solves det(P-V_0, V_n-V_0, ..., V_1-V_0)=0, i.e. every point is on hyperplane
1126     see Bower & Woodworth, A programmer's geometry, Butterworths 1983.
1127 
1128   derivation of 3-d minnorm
1129     Goal: all vertices V_i within qh.one_merge of hyperplane
1130     Plan: exactly translate the facet so that V_0 is the origin
1131           exactly rotate the facet so that V_1 is on the x-axis and y_2=0.
1132           exactly rotate the effective perturbation to only effect n_0
1133              this introduces a factor of sqrt(3)
1134     n_0 = ((y_2-y_0)*(z_1-z_0) - (z_2-z_0)*(y_1-y_0)) / norm
1135     Let M_d be the max coordinate difference
1136     Let M_a be the greater of M_d and the max abs. coordinate
1137     Let u be machine roundoff and distround be max error for distance computation
1138     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
1139     The max error for distance of V_1 is sqrt(3) u M_a M_d M_d / norm.  Offset=0 at origin
1140     Then minnorm = 1.8 u M_a M_d M_d / qh.ONEmerge
1141     Note that qh.one_merge is approx. 45.5 u M_a and norm is usually about M_d M_d
1142 
1143   derivation of 4-d minnorm
1144     same as above except rotate the facet so that V_1 on x-axis and w_2, y_3, w_3=0
1145      [if two vertices fixed on x-axis, can rotate the other two in yzw.]
1146     n_0 = det3_(...) = y_2 det2_(z_1, w_1, z_3, w_3) = - y_2 w_1 z_3
1147      [all other terms contain at least two factors nearly zero.]
1148     The max error for n_0 is sqrt(4) u M_a M_d M_d / norm
1149     Then minnorm = 2 u M_a M_d M_d M_d / qh.ONEmerge
1150     Note that qh.one_merge is approx. 82 u M_a and norm is usually about M_d M_d M_d
1151 */
qh_sethyperplane_det(qhT * qh,int dim,coordT ** rows,coordT * point0,boolT toporient,coordT * normal,realT * offset,boolT * nearzero)1152 void qh_sethyperplane_det(qhT *qh, int dim, coordT **rows, coordT *point0,
1153           boolT toporient, coordT *normal, realT *offset, boolT *nearzero) {
1154   realT maxround, dist;
1155   int i;
1156   pointT *point;
1157 
1158 
1159   if (dim == 2) {
1160     normal[0]= dY(1,0);
1161     normal[1]= dX(0,1);
1162     qh_normalize2(qh, normal, dim, toporient, NULL, NULL);
1163     *offset= -(point0[0]*normal[0]+point0[1]*normal[1]);
1164     *nearzero= False;  /* since nearzero norm => incident points */
1165   }else if (dim == 3) {
1166     normal[0]= det2_(dY(2,0), dZ(2,0),
1167                      dY(1,0), dZ(1,0));
1168     normal[1]= det2_(dX(1,0), dZ(1,0),
1169                      dX(2,0), dZ(2,0));
1170     normal[2]= det2_(dX(2,0), dY(2,0),
1171                      dX(1,0), dY(1,0));
1172     qh_normalize2(qh, normal, dim, toporient, NULL, NULL);
1173     *offset= -(point0[0]*normal[0] + point0[1]*normal[1]
1174                + point0[2]*normal[2]);
1175     maxround= qh->DISTround;
1176     for (i=dim; i--; ) {
1177       point= rows[i];
1178       if (point != point0) {
1179         dist= *offset + (point[0]*normal[0] + point[1]*normal[1]
1180                + point[2]*normal[2]);
1181         if (dist > maxround || dist < -maxround) {
1182           *nearzero= True;
1183           break;
1184         }
1185       }
1186     }
1187   }else if (dim == 4) {
1188     normal[0]= - det3_(dY(2,0), dZ(2,0), dW(2,0),
1189                         dY(1,0), dZ(1,0), dW(1,0),
1190                         dY(3,0), dZ(3,0), dW(3,0));
1191     normal[1]=   det3_(dX(2,0), dZ(2,0), dW(2,0),
1192                         dX(1,0), dZ(1,0), dW(1,0),
1193                         dX(3,0), dZ(3,0), dW(3,0));
1194     normal[2]= - det3_(dX(2,0), dY(2,0), dW(2,0),
1195                         dX(1,0), dY(1,0), dW(1,0),
1196                         dX(3,0), dY(3,0), dW(3,0));
1197     normal[3]=   det3_(dX(2,0), dY(2,0), dZ(2,0),
1198                         dX(1,0), dY(1,0), dZ(1,0),
1199                         dX(3,0), dY(3,0), dZ(3,0));
1200     qh_normalize2(qh, normal, dim, toporient, NULL, NULL);
1201     *offset= -(point0[0]*normal[0] + point0[1]*normal[1]
1202                + point0[2]*normal[2] + point0[3]*normal[3]);
1203     maxround= qh->DISTround;
1204     for (i=dim; i--; ) {
1205       point= rows[i];
1206       if (point != point0) {
1207         dist= *offset + (point[0]*normal[0] + point[1]*normal[1]
1208                + point[2]*normal[2] + point[3]*normal[3]);
1209         if (dist > maxround || dist < -maxround) {
1210           *nearzero= True;
1211           break;
1212         }
1213       }
1214     }
1215   }
1216   if (*nearzero) {
1217     zzinc_(Zminnorm);
1218     /* qh_joggle_restart not needed, will call qh_sethyperplane_gauss instead */
1219     trace0((qh, qh->ferr, 3, "qh_sethyperplane_det: degenerate norm during p%d, use qh_sethyperplane_gauss instead.\n", qh->furthest_id));
1220   }
1221 } /* sethyperplane_det */
1222 
1223 
1224 /*-<a                             href="qh-geom_r.htm#TOC"
1225   >-------------------------------</a><a name="sethyperplane_gauss">-</a>
1226 
1227   qh_sethyperplane_gauss(qh, dim, rows, point0, toporient, normal, offset, nearzero )
1228     given(dim-1) X dim array of rows[i]= V_{i+1} - V_0 (point0)
1229     set normalized hyperplane equation from oriented simplex
1230 
1231   returns:
1232     normal (normalized)
1233     offset (places point0 on the hyperplane)
1234 
1235   notes:
1236     if nearzero
1237       orientation may be incorrect because of incorrect sign flips in gausselim
1238     solves [V_n-V_0,...,V_1-V_0, 0 .. 0 1] * N == [0 .. 0 1]
1239         or [V_n-V_0,...,V_1-V_0, 0 .. 0 1] * N == [0]
1240     i.e., N is normal to the hyperplane, and the unnormalized
1241         distance to [0 .. 1] is either 1 or   0
1242 
1243   design:
1244     perform gaussian elimination
1245     flip sign for negative values
1246     perform back substitution
1247     normalize result
1248     compute offset
1249 */
qh_sethyperplane_gauss(qhT * qh,int dim,coordT ** rows,pointT * point0,boolT toporient,coordT * normal,coordT * offset,boolT * nearzero)1250 void qh_sethyperplane_gauss(qhT *qh, int dim, coordT **rows, pointT *point0,
1251                 boolT toporient, coordT *normal, coordT *offset, boolT *nearzero) {
1252   coordT *pointcoord, *normalcoef;
1253   int k;
1254   boolT sign= toporient, nearzero2= False;
1255 
1256   qh_gausselim(qh, rows, dim-1, dim, &sign, nearzero);
1257   for (k=dim-1; k--; ) {
1258     if ((rows[k])[k] < 0)
1259       sign ^= 1;
1260   }
1261   if (*nearzero) {
1262     zzinc_(Znearlysingular);
1263     /* qh_joggle_restart ignored for Znearlysingular, normal part of qh_sethyperplane_gauss */
1264     trace0((qh, qh->ferr, 4, "qh_sethyperplane_gauss: nearly singular or axis parallel hyperplane during p%d.\n", qh->furthest_id));
1265     qh_backnormal(qh, rows, dim-1, dim, sign, normal, &nearzero2);
1266   }else {
1267     qh_backnormal(qh, rows, dim-1, dim, sign, normal, &nearzero2);
1268     if (nearzero2) {
1269       zzinc_(Znearlysingular);
1270       trace0((qh, qh->ferr, 5, "qh_sethyperplane_gauss: singular or axis parallel hyperplane at normalization during p%d.\n", qh->furthest_id));
1271     }
1272   }
1273   if (nearzero2)
1274     *nearzero= True;
1275   qh_normalize2(qh, normal, dim, True, NULL, NULL);
1276   pointcoord= point0;
1277   normalcoef= normal;
1278   *offset= -(*pointcoord++ * *normalcoef++);
1279   for (k=dim-1; k--; )
1280     *offset -= *pointcoord++ * *normalcoef++;
1281 } /* sethyperplane_gauss */
1282 
1283 
1284 
1285